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SUMMARY 

A computer program has been developed which will calculate the particle 
spectra and dose behind multilayer, infinite slabs of shielding from any 
given energy spectrum of protons impinging normally on the shield. Both 
spectra and dose are calculated for the incident primary protons that pene- 
trate the shield and for the following types of secondary radiation produced 
in the shield: cascade protons, cascade neutrons, and evaporation neutrons. 

This program is written in Fortran IV for the IBM 7094 II computer. 


INTRODUCTION 

Space vehicles may require shielding to protect astronauts and radiation 
sensitive components from the hazards of space radiation, the most hazardous 
of which appears to be the proton radiation. The important sources of proton 
radiation are the Van Allen belts, solar flares, and galactic cosmic radiation. 
When space vehicles are irradiated by protons the interaction of protons with 
atomic nuclei in the vehicle produces secondary nucleons. The most important 
of these secondary nucleons appears to be the cascade protons, cascade neutrons, 
and evaporation neutrons. In order to evaluate the shielding requirements it 
is necessary to evaluate the dose caused by primary and secondary radiation. 

Previous computer programs on space radiation shielding have been written 
by (at least) three groups (see refs. 1 to 3). 

References 1 and 2 used the 1958 data of Metropolis et al. for estimating 
the dose from secondary radiation. The reference 3 program was not designed 
for general distribution. Recently H. Bertini of ORNL (ref. 4) developed a 
Monte Carlo program (with an improved nuclear model compared to previous work) 
for estimating the yields of secondaries and their energies. These reasons 
combined with the fact that the Lewis computer is now restricted to Fortran TV 
provided the motivation at the Lewis Research Center to develop a computer 
program capable of calculating the doses from primary proton and secondary 
radiation produced in the shield. 

The Lewis program (called LPSC) calculates the particle spectra and dose 
inside a multilayer slab shield due to protons impinging normally on the outer 
face of the shield. The slabs are infinite in extent and have finite thick- 
ness. Spectra and doses may also be calculated at intermediate thicknesses 
through the shield known as print bounds. 
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The total dose evaluated includes the doses from primary protons and 
the secondary radiation consisting of cascade protons, cascade neutrons, and 
evaporation neutrons. Evaporation protons contribute a negligible amount 
to the total dose. Other secondary radiation (heavy particles, pions, 
secondary gammas, etc.) are not evaluated by this program. The flux spectra 
of the primary protons, cascade protons, and cascade neutrons are calculated 
in this program. 

This program calculates the penetration of normally incident protons, 
cascade protons, and cascade neutrons using a straight ahead approximation 
(the primary particle and cascade secondaries are emitted in the same di- 
rection as the primary incident particle) . This approximation has been 
shown by reference 5 to be valid for protons with energies >1D0 MeV. For 
energies < 100 MeV this method may over estimate the dose. The evaporation 
neutrons were assumed to be emitted isotropically. 

It may be shown, using the straight-ahead approximation, that the dose 
received at the center of a sphere due to an isotropic flux outside of the 
sphere is the same as the dose received behind a slab shield of the same 
thickness where a flux of the same magnitude is impinging normally to the 
slab. This will also be true for doses from those secondary particles as- 
sumed to be emitted in the direction of the primary particles. Thus the 
slab doses calculated for the normally incident primary protons and the 
cascade secondaries will be the same as the doses at the center of the 
sphere for an isotropic flux outside; however, the evaporation neutron doses 
are not simply related and pertain only to the slab geometry. 

A description of the program and program listing are found on pages 
16-32 ••and 46r-73. 


SYMBOLS 


Some of the symbols used in this report are listed below. Other symbol 
are defined in the text. 


D 


cn 


D, 


cp 


D, 


en 


D. 


PP 


dose from cascade neutrons, rad or rad/hr and rem or rem/hr 
dose from cascade protons, rad or rad/hr and rem or rem/hr 
dose from evaporation neutrons, rad or rad/hr and rem or rem/hr 
dose from primary protons, rad or rad/hr and rem or rem/hr 
energy in MeV 


G(X k ,9 m ) attenuation kernel for evaporation neutrons see text for dimen- 
sions 
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N(>E) 


E(>P) 

N(E ± ,X) 

V'j'V 

N cn (Ej,X 4 ) 

c M p (I., E;) ,8X k ) 


e H p (Ei,SX k ) 

A( E i> E J.8X k ) 

e Bn(Ei.8X k ) 

= P p (S 1 ,E.,8X k ) 

’8X k ) 


X 


y cpp( E i' E j) 


number of protons with energy greater than E, protons /cm 2 
or protons /cm 2 sec 

number of protons with rigidity greater than P, 
protons/cm 2 or protons/cm 2 sec 

the primary proton flux at energy E^ at depth X, 
protons/cm 2 or protons/cm 2 sec 

cascade proton flux with energy Ej which penetrates the 
shield, protons/cm 2 or protons/cm 2 sec 

cascade neutron flux with energy Ej which penetrates the 
shield, neutrons/cm 2 or neutrons/cm 2 sec 

cascade neutrons with energy Ej produced by incident pro- 
tons with energy Ep in layer 5X k , neutrons /cm 2 or 
neutrons / cm 2 sec 

evaporation neutrons produced by protons with energy Ep 
in layer 5X k , neutrons/cm 2 or neutrons/cm 2 sec 

cascade neutrons with energy Ej produced by neutrons with 
energy Ep_ in layer &X k , neutrons/cm 2 or neutrons/cm 2 sec 

evaporation neutrons produced by neutrons with energy Ep 
in layer 6X k , neutrons /cm 2 or neutrons/cm 2 sec 

cascade protons with energy Ej produced by protons with 
energy Ep in layer 6X k , protons/cm 2 or protons/cm 2 sec 

cascade protons with energy Ej produced by neutrons with 
energy Ep in layer 5X k , protons /cm 2 or protons/cm 2 sec 

depth into the shield, g/cm 2 

yield of cascade protons with energy Ej per interacting 
proton with energy Ep 

yield of cascade neutrons with energy Ej per interacting 
proton with energy Ep 

yield of evaporation neutrons per interacting proton with 
energy Ep 


yenp( E i) 
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y (e.,eA 

•'cpn i’ y 


^enn^ E i^ 

n(%,X) 


yield of cascade protons with energy Ej per interacting 
neutron with energy 

yield of cascade neutrons with energy Ej per interacting 
neutron with energy E^ 

yield of evaporation neutrons per interacting neutron with 
energy Ej_ 

neutron flux with energy E^ at depth X, neutrons/cm^ or 
neutr ons/ cm? sec 


rigidity, MV 


METHOD OF CALCULATION 
Introduction to Method of Calculation 

The initial flow of a calculation progresses as follows; see figure 1. 
The computer reads in a maximum incident energy E(Max) in_MeV at (1), Using 
the proton range energy data an E(Max) at the exit face ( g)' is calculated. 
The desired energy group bounds E]_, Eg, • . . E^, E(Max) are read in and 
assumed at the exit face (3). AE values are calculated by AEp =* (E^ - 0), 
AEg •* (Eg - Ep), . . ., AEn » (E(Max) - En) » The number of small 8E's con- 
tained in each AEp are read. The AE^’s and BE's enable the computer to 
calculate the energy bounds on the . exit face. Then the energy bounds sire 
calculated on the incident face (' 4 ) . The number of protons in each energy 
group are calculated and the production of secondary particles and the at- 
tenuation of all particles for each group is accomplished. 

Figure 1 shows three print bounds at B, C, and D where data are printed. 
Incident proton dose is calculated at the incident face A also. 


Proton Energy Groups and Spectra 

This program allows the user a choice of reading in the incident energy 
boundaries or the exit energy boundaries. The use of exit energy boundaries 
is to be preferred for running a proton energy spectrum. Reading in the 
incident energy boundaries are preferred for a monoenergetic case. 

Some of the details for the case where exit energy boundaries are assumed 
are as follows: 
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Tables of range versus energy are supplied as input data for the shield 
materials being investigated. The Lagrange 2 point interpolation equation 
is used to interpolate in the tables. The program selects the maximum 
incident energy E(Max) and calculates the range associated with this energy 
for the first shield material encountered, R(E(MAX), Material l) = R(E,l). 

The thickness T(l) of this first layer is then substracted from R(E,l). 

The residual range r(E(l),l) of this proton group after passing through 
T(l) is given by r(E(l),l) « R(E,l) - T(l). The energy E(l) associated 
with r(E(l),l) is calculated and the range for a proton of this energy in 
the next material is calculated. The thickness of this layer is then sub- 
tracted off and this process is repeated until the maximum energy at the 
exit face of the shield has been computed. 

When the maximum energy at the exit face of the shield is known, the 
program calculates the energy bounds at the exit face according to the input 
data. There are two degrees of subdivision for the energy bounds. AE 
represents a large interval and BE represents a small interval. The AE 
boundaries axe used as convenient points at which the magnitude 6E(j) is to 
change . 

The interval bounds of AE and the number n of 6E(«J) in each AE 
at the exit face are input data. The 8E(j) values are calculated by the 
program. The energy bounds E(l) and energy intervals 6E(l) on the incident 
face corresponding to the energy bounds E(j) and energy intervals 6E(j) 
on the exit face are then calculated. The average energy of each group I 
is calculated by. 


I(l) . E(I+ ., 1) + , S(I) (I) 

2 

The average energy E(l) is used to calculate the cross sections, 
yields, and energy of secondary paxticles. The average energy incident 
on the next lajrer is obtained by degrading* thfe'^-energy bounds^ an<l' palcu* • 
lating a new E(l). This is due to the nonlinear change in BE which 
makes E(l) degraded not equal to the new average energy. 

The incident energy associated with zero energy at the exit face of the 
shield is not the minimum incident energy considered because lower energy 
protons, although they do not penetrate the shield, may produce secondary 
neutron radiation that can penetrate the shield. The procedure used for 
finding incident energy intervals below this energy is illustrated by the 
following example. 

Figure 1 shows four print bounds A, B, C, and D. At D, E(0) « 0 MeV 
is the zero energy point bn the exit face of the shield. Similarly there 
exists a zero energy point at each print bound. At print bound C the 
energy interval from E(0) « 0 to C subdivided by the program in 
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one of two ways, (l) This interval may be subdivided the same as this 
corresponding interval on the exit face D or a new set of subdivisions 
may be read in. The corresponding intervals on the incident face are 
found. (2) This interval may be divided into an equal number of increments. 
The other print bounds wijl be handled in the same manner . 

It was found by running the program that if the same energy bounds were 
used for neutrons as for protons the neutrons tended to cluster into a few 
energy groups. Therefore to improve the distribution of dose with energy 
separate energy bounds were required for neutrons and protons. The neutron 
energy bounds are read in as a table. 

The incident proton spectra may be simulated by one of the following 
equations. Let either N(>E) or dN/dE be given by r(E) as indicated in 
equations (2) and (3). 

r(E) « N(>E) (2) 

or 

T(E) ■ dE/dE (3) 

then r(E) , the proton spectrum, may be calculated by any one of equations 
(4), (5), (6), or (7). 


r(E) -A E' b 

T(E) = A(E)exp(-B(E)) 

4 

Log (r(E)) - A* E 1-1 

i-1 


Log (r(E)) 


2 


Aj_(Log E) 1-1 


(4) 

(5) 

( 6 ) 

(7) 


Where the constants or Aj. are selected compatible with either 

equation (2) or (3). See page 26 for a definition of ' A'(E) ahd 
B(E). The program will also accept rigidity spectra where 


then 


r(p) = N(>P) 

r(P) * dN/dP 


( 8 ) 

(9) 
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r(P) = A exp (-P/P Q ) (10) 

If equations (2) through (10) fail to represent the desired spectrum a 
table of values may be read in for N(>E) or dN/dE as a function of energy. 

When the spectrum is read in as a table the values of dN/dE or 
N(>EjJ, whichever the table contains, will be read and/or interpolated for 
each value of E-^ required by the program. 

When the choice of spectrum (equation or table) has been made, the 
continuous input spectrum is approximated by selecting a finite number of 
proton energy groups. The number of protons in each group is calculated. 

For the differential spectrum the number of protons in the i th group is given 

by 


N(E ,0) 



( 11 ) 


The width of the i th energy interval is SEj.. The midpoint of the i th inter- 
val is Ej. The quantity (dN/dE)_ is calculated from one of the equations 
or a table. ®i 

For the integral spectrum, the number of protons in the i th group is 
given by 


N(E 1 ,0) ~ N(>E i ) - N(>E 1+1 ) (12) 

When running a monoenergetic case one should select a finite interval 
width 8E such that the midpoint of 5E is the required energy %p. A 
small 8E on the order of 10“ 2 MeV is recommended here. The calculation 
is not very sensitive to 5E where 10”3 < 5E < 0.1 MeV. For 5E < 10~3 
significant figures are lost during substraction. 


Proton Attenuation and Production of Secondaries 

When protons penetrate a shield energy is lost due to ionization inter- 
actions. The range in a given material and rate of energy loss per unit path 
length are energy dependent. Protons that are stopped by ionization are re- 
moved from the beam. Incident protons which experience inelastic collisions 
with nuclei in the shield are removed from the beam. The secondary protons 
and neutrons produced in these inelastic collisions are added to the beam. 

Secondary particles are produced by primary protons, cascade protons, 
and cascade neutrons. Secondaries produced by evaporation neutrons are not 
calculated by LPSC. 
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The secondary particles are calculated as follows: The various slabs of 

materials in the shield are divided into increments 8X k thick and print 
boundaries Xq (see fig. 2). The print bounds are locations within the shield 
at which data are printed. The SXfc are used as secondary source regions 
throughout the shield. Between two consecutive print bounds all SX k values 
are the same. When a print bound is crossed the 6X k are again all equal 
but they may be different (in number and/or size) from those in the previous 
region. 

Indexing presents a problem in a calculation of this nature. In general 
the index i or I is used to indicate incident particle energies and j 
or J for exit or secondary particle energies. At an internal boundary the 
energies incident on the K^ 1 layer are the same as those that exit the K-l 
boundary. However, even though these energies are the same when viewed as 
incident particles the subscript i or I will be used and when viewed as 
exit particles the subscript j or J will be used. One exception will be 
found in the dose calculation where the exit primary proton energies retain 
the index i to prevent confusion with secondary particle energies which use 
the subscript j. These primary particles that exit the shield may be con- 
sidered incident on a detector to preserve the above rules. 

Let N(E^,X) be the number of protons / cm^ (or protons/ cm 2 sec) in each 
energy group Ej_ at depth X into the shield which are incident on a layer 
6X k> Let N(Ej,X + 8X k ) represent the number of protons at energy Ej from 

N(E^,X) which pass through SK^ without experiencing a nuclear interaction. 

A drop in energy from E^ to Ej occurs due to ionization. Let Lp(Ei) be 
the macroscopic cross section in cm^/g for proton inelastic collisions. The 
two quantities N(Ei,X) and N(Ej,X + SX) are related by 

N(Ej,X + 8X) « N(E ± ,X) exp (-EpC^) * 5X k ) (l3) 

Since exp (-Sp(E^) ‘ 8X k ) is the probability that no interactions occur in 
5X k then [l - exp(-£p(E^) »&X k )] is the probability that an interaction does 
occur in 8X k . Hence the number in E(E^,X) that interact in BX k is given 
by N(Ei,X)[l - exp(-Ep(E)5X k ) ] . Let y cpp ( E i^ E j) be the average yield of 
cascade protons at energy Ej produced by an interacting proton at energy 
E^. Let c Pp(E^,Ej,5X k ) be the cascade protons at energy Ej produced by 
protons of energy Ej in layer 8X k . The secondary group is related to 
the incident group by 

c P A'V“k> - - e*p[-£ p (E.)8X k )} x cpp( S i’V 


( 14 ) 
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Similarly let c Ep(E^,Ej ,SX k ) be the cascade neutrons and ycnp(%>®j) 

the average yield of cascade neutrons. Then the cascade neutrons. are re- 
lated to the incident protons by 

cN p (Ei, Ej ,5X k ) - Nd-pX) jl - exp[-Ep(%)5X k ]j y cnp (I;pE.) 

(15) 

Let e N p (E.£,8X k ) be the yield of evaporation neutrons produced by 
the incident proton group N(!j^,#)* Let yenp(%) be ^he average yield of 
evaporation neutrons per interacting proton at energy E^. The evaporation 
neutrons are assumed to have a fission spectrum. The evaporation neutrons 
produced are related to the incident protons by 

e N p (Ep6X k ) « N(EpX) |l - exp[-S p (K l )6X k ]j y^E^ (16) 

When neutrons are incident on a layer SX k with sufficient energy E^ 
to produce secondaries a similar set of equations can be written 

c P J E l’ I V SX k> - "< 5 i- X > f 1 - ycpn(%> E j) 

cN n (ii,Ej,6X k ) - nd^X) [l - exp[-E n (S 1 )Sx k ]] y^di.Ej) 

(18) 

e N n<V 6X k> - n < 5 i> X > f 1 - expI-E^)^]] y^ffp (19) 

where n(Ej^,X) represents the incident neutrons/cm^ (or neutrons/cm^ sec) at 
energy Ej_ at depth X and E n (Ej_) represents the neutron inelastic cross 
sections in cm^/g. Equations (17), (18) and (19) give the cascade protons, 
cascade neutrons, and evaporation neutrons, respectively. 

The secondary particles produced in a layer SX k are placed at the 

center of this layer for calculating attenuation through the remainder of 
the shield. These secondaries are attenuated across the half layer 5X k /2 
(in which they were born) by removing those that experience inelastic 
collisions. No higher generation secondary production is calculated in this 
half layer due to the secondaries that were born here. 

Particle penetration through each 5X k is accomplished using the 
straight ahead approximation except for evaporation neutrons which are as- 
sumed to be emitted isotropically. 
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The data from reference 4^ gives secondary yields, energies of secondaries, 
and inelastic cross sections at energies . 25 MeV < E < 400 MeV for C, 0, Al, 

W, Pb, and U for protons and neutrons bombarding. Data for N, Ti, and Fe 
were obtained by interpolating the data from reference 4 as a function of mass 
number. 

The secondary particle data tables in LPSC contain data up to 1000 MeV. 

The secondary yields and energy of secondaries were extrapolated from 400 MeV 
to 1000 MeV by fairing in a curve. For neutron and proton energies > 400 MeV 
the inelastic cross sections at 400 MeV are used. Low energy (<25 MeV) neutron 
cross section data were taken from the literature. 

The proton cross section tables contain entries in the energy range 
10 MeV to 1000 MeV for hydrogenous materials. For nonhydrogenous materials 
the proton cross section tables contain data in the energy range 25 MeV to 
400 MeV. When data beyond the range of the tables are required for yields 
and energy of secondaries the program will extrapolate the tables using a 
2 point Lagrange interpolation equation. If neutron cross sections are called 
for below the minimum value in the table the cross section is set equal to 
zero. The proton range energy table contains data up to 10® MeV. 

The yields for cascade neutrons emitted when cascade neutrons are in- 
cident were carried to the low energy threshold for inelastic scattering by 
assuming the yield below 25 MeV would be given by the ratio cr n n '/o u 
Where CT n n ' is the inelastic scattering cross section and cf n x is the total 
inelastic cross section. This was accomplished for the following materials in 
this program: carbon, oxygen, nitrogen, aluminum, titanium, iron, uranium, 

water and polyethylene. For lead and tungsten the n, 2n cross section was 
included for energies < 15 MeV. 


Attenuation of Secondaries 

The cascade protons are attenuated in the same manner as the primary 
protons . 

The straight ahead approximation is used for both cascade protons and 
cascade neutrons. Additional generations of secondaries produced by second- 
aries for both cascade protons and cascade neutrons are calculated or deleted 
on command. 

In this program cascade neutrons are assumed to experience inelastic 
collisions with nuclei having mass numbers >9. This interaction assumes that 
the incident neutron was absorbed followed by the emission of secondary cascade 


1 


Reference 4 also gives data for other elements not used in LPSC. 
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and evaporation particles. This is the same type interaction described for 
protons only now the bombarding particles are neutrons. Elastic collisions 
of cascade neutrons with nuclei heavier than the neutron tend to produce 
small deflections which result in small energy loss, hence a small attenu- 
ation. 

Elastic collisions of cascade neutrons in hydrogen do have a significant 
effect on attenuation. This results because the mass of the target is nearly 
equal to the mass of the incident particle. Some simplifying assumptions 
were made regarding the elastic collision process in hydrogen. The energy of 
a neutron scattered off a hydrogen nucleous was averaged (and weighted using 
differential cross sections) over all angles of scattering. The recoil 
nucleous was assigned an energy which was the difference between the incident 
neutron energy and the scattered neutron energy. Both the neutron and proton 
were assumed to be emitted in the direction of the incident neutron. Tables 
of scattered neutron energies and recoil nucleous energies were calculated 
for several incident neutron energies. 

The elastic collision of a proton on hydrogen nuclei was assumed to 
yield two protons both having the same direction as the incident particle 
and each proton having half the energy of the incident particle., 

The above assumptions for neutrons and protons enables data tables to 
be constructed for hydrogen similar to the secondary_ particle data tables 
for the other materials. 

The method used to attenuate the dose from evaporation neutrons is 
similar to the one described in reference 2. The details of this calcula- 
tion are presented in the section on dose calculations. 

The data tape for the LPSC program contains beryllium as the only choice 
of shield material which does not contain secondary yields. The incident- 
proton spectrum is attenuated by ionization only. 


Primary and Cascade Proton Dose 


The primary proton flux N(E^,X|) which penetrates the shield is con- 
verted to dose in rads by multiplying by the stopping power dE/dX (for water 
or tissue) and a unit conversion factor U. U _is in rad g/MeV or 
rad g sec/MeV hr depending on the units of N(E^,X|). Let the primary pro- 
ton rad dose be given by Dpp then 


D. 


PP 


= U 





( 20 ) 
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The rem dose is given by 

I 

Spp - 0 ^ N < f i> x ^(S) ii RBE '%> < 21 > 

i=l 

The RBE factors as a function of proton energy are inserted in the program 
as a table The RBE data was obtained from reference 6. 

Equations (20) and (2l) are used to calculate the dose in rad and_rem, 
respectively, for cascade protons where E(E i ,Xg) is replaced by N C p(Ej,X^), 

summation is accomplished over the index j which represents cascade proton 
energies, and ( dE/dX) „ . 

Since dE/dX increases with a decrease in E there exists a possibility 
of a given proton spectrum producing a higher dose at some positive depth than 
existed at the incident face. In order to observe this effect the increase in 
dE/dX must outweigh the loss in particle attenuation. 


Cascade Neutron Dose 

_If the cascade neutron flux that penetrates the shield is given by 

N (E.,X„) and the cascade neutron dose is given by D then 

cn j £ e J cn 


D, 


cn 



N cn(Ej,Xg)A(Ej) 


( 22 ) 


J=1 


where A(Ej) is the flux to dose conversion factor at neutron energy Ej . 

The values of A(eO for the energy interval 0.1 MeV to 10 MeV were ob- 

d 

tained from reference 6. The values of A(E.) for the energy interval 60 MeV 

J 

to 400 MeV were obtained from reference 7. The data from 10 to 60 MeV were 
faired in. These conversion factors are based on calculated maximum doses 
produced in a slab of tissue by normally incident neutrons. 


Evaporation Neutron Dose 

The evaporation neutron dose is calculated in a manner similar to that 
of reference 2 with some alterations. The dose, due to evaporation neutrons 
born in each SX^, is evaluated as the dose from an infinite plane source at 
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the center of 8X^ (see fig. 3). The Albert Welton Kernel is integrated over 
angles and shield layers which are neutron sources. 

The rad dose from evaporation neutrons in the entire slab is given by 


D en = 



+ e N n( E ^ 5X k) 

J-l 




(23) 


where 


8fl ( 0 m^ 0 m+l) = 2lt ( cos 0 m - cos 0 m+1 ) 


(24) 


The index m represents the number of angles measured in the interval 
0 < 0^ < 90°, see figure 3. The angle 6^ is the midpoint of the angular 
interval (® m #® m +i) • The index k represents the number of increments of 
shield thickness SX^. The index i represents the number of energy groups 

of incident protons which produce evaporation neutrons in 6X^. The index j 
represents the number of energy groups of incident cascade neutrons which pro- 
duce evaporation neutrons in SX-^. The units for these sums are particles/cm 2 

or particles/cm 2 sec. 

The maximum number of angles that can be used is 10. Experience has 
shown that five angles areL adequate' in reproducing' the dose caletdationsito 
two significant figures when compared with the 10 angle calculation. For 
materials where the evaporation neutron dose is small it may be preferred to 
run with one angle because the running time is less. 

The term G(X^^0^ 1 ) is the attenuation kernel. For hydrogenous material 
G is the Albert Welton kernel using the coefficients derived by Casper 
(ref. 8). For nonhydrogenous material G is used as shown by reference 2. 

The function G(X k ,0 1 { 1 ) ’ is calculated in LPSC as follows 



( 25 ) 


14 


where 


F(t)) = T] exp( - C 3 T] ), TJ > 2.0 


(26) 


Equation (27) represents a straight line extrapolation of F(t)) versus T) for 
0 < T) < 2.0. 


F(tj) » 0.772 - 0.065 tj, 0 < T] < 2.0 
^max 


I'i 


n=k 


S n = removal cross section for material n in cm^/g 
r n = slant path length through material n in g/crn^ 


K n = a constant 


( 27 ) 


( 28 ) 


K n = 1.0 for all hydrogenous materials 

K n - 1.0 for 2 < Z < 6 for nonhydrogenous materials 

K n = 0.5 for Z > 6 for nonhydrogenous materials 

H n = the ratio of hydrogen density in material n to hydrogen density in 
water 

P n = the density of material n in g/cm^ 

Kmax ~ the maximum value to the index on the number of SX n layers being 
calculated to the detector point. As the detector point moves the 
value of K max will change. 

The index in equation (25) starts at K and progresses to K^y . This 
indicates that source layers are numbered 1 starting at the incident face and 
progress to K^y at the detector face. Therefore the attenuation applied to 
the layer is from K to K^y . As the calculation progresses to a new 
print bound the value of will change. See figure 2. 

(C 1 ,C 2 ,C 3 ,C 4 ) = (5.389 X 10 -9 , 0.3492, 0.4223, 0.6984) (29) 

Casper in reference 8 derived the coefficients in equation (29) to fit the 
data in the range of 10 cm to 130 cm from the source in water. 
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■ Equation (29) represents a departure from reference 2. 

The value of = 5.389 X.' 10 “^ is for source terms in particles/cm^. 

Then G is in ( rad/flare)/( neutrons /cn£ ) . If the source terms are in 
neutrons/cm 2 sec then Cp = 3600 X 5.389 x: 10 -9 and G is in (rad/hr)/ 
(neutrons/cm^ sec). 

If the shield is all nonhydrogenous G^^.,©^) is calculated using 
equations (25) and (27) where T] = 0. If the shield is all hydrogenous 
G(Xj t ,0^) is calculated using equations (25), (26), and (27) where rj ^ 0. 

When some layers following layer n contain hydrogen and some do not, 
the method of reference 2 is used. Replace S n with S n Z n . Then if non- 
hydrogenous material follows layer n, set Z n = K n . If hydrogenous material 
follows layer n, then Z n is selected as follows: 

If 

K ■ 1.0 then Z n = 1.0 (30) 

if 

K n = 0.5 then Z n is taken as the minimum of 


Z n = 1.0 or 
K 

'max 



n=k+l 


The Albert-Welton kernel (eq. 25) with F(rj) defined as in equation (26), 
pertains to neutrons with a fission energy spectrum. Bertini's data indicates 
that the evaporation neutron energy spectrum is harder than the fission spectrum 
when incident particle energies are > 25 MeV. This will tend to make the calcu- 
lated dose from these evaporation neutrons low. Evaporation neutrons do not 
produce secondary particles in this program. 

Equation (23) calculates the rad dose from evaporation neutrons. The rem 
dose is obtained by multiplying the rad dose by an RBE of 10. The RBE can be 
readily changed in the program to any value considered more applicable. 


Total Doses and Flux and Source Terms 

In addition to calculating the individual dose components previously 
mentioned the program calculates the total proton dose, the total neutron dose. 
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and the total dose from all particles. If dN/dE or N(>E) is input in 
protons/(cm^ MeV flare) or protons/(cm2 flare), the doses calculated are 
rad/flare and rem/flare. If dN/dE is in protons /(cm? MeV sec) then the 
doses calculated are in rad/hr and rem/hr. 


The LPSC will calculate the particle flux at each print bound in 
particles/cm^ or particles/cm^ sec. The evaporation neutron source strengths 
in units of neutrons/g or neutrons/g sec are also calculated for each SX 
layer . 


PROGRAM DESCRIPTION LPSC 

The Lewis Proton Shielding Code (LPSC) consists of the main program 

PISR and the following subroutines: 

FLUXEQ calculates initial proton flux as a function of initial incident 
energy. 

INVALU sets up energy intervals at exit face of slab and at all inter- 
mediate print out bounds at which printed output is desired; 
uses these energy groups to calculate initial incident energy 
intervals and with FLUXEQ establishes initial proton spectrum. 

EVNEDO calculates the evaporation neutron dose based on the source 
terms from the midpoints of each SX increment. 

XS computes the cross sections of protons and neutrons as a function 

of energy. 

LAGRNG interpolation scheme based on Lagrange fundamental formula for 
interpolation . 

YIELDS calculates the yield of secondary particles per collision as a 
function of the type and energy of the bombarding particle. 

RANGE a dual purpose subroutine which calculates the range as a function 
of energy, and the energy as a function of range. 

CASNRG computes the energy of cascade protons and neutrons as a function 
of the type and energy of the bombarding particle. 

DOSEK computes proton and neutron flux to dose conversion factors for 
doses in rad (or rad/hr) and rem (or rem/hr) as a function of 
energy. 

PROPTY transmits all material properties from magnetic tape to disk 

storage for faster access and transmits tables of flux to dose 
conversion factors to core storage. 


17 


SORT is a general purpose sorting routine. 

The main program PISR is divided into four separate sections. In the 
first section the input data is entered into the computing machine and large 
blocks of storage are initialized for later use. The second section contains 
the calculations for the attenuation of the primary protons, and the colli- 
sions producing secondary particles. All secondary particles are assumed to 
be born at the midpoint of each SX subinterval and are attenuated across 
the second half of the 8X sub interval. The dose calculations are contained 
in the third section of the main program, and section four controls the output 
of data. 

Subroutine FLUXEQ calculates the flux values of the proton spectrum, 
incident on the shield at various energies. In subroutine INVAHJ an equa- 
tion number, entered as input, specifies the particular analytical equation 
that simulates the spectrum or the table of values of flux versus energy that 
defines the spectrum. These equations represent either the integral or dif- 
ferential spectral forms . The .code converts the rigidity spectral equations 
to energy and all calculated spectra are presented in terms of energy. 

The INVAHJ subroutine is a two section program which computes all the 
initial data for a particular problem. The first section contains the com- 
putation of the AX intervals and 5X subintervals and then establishes 
the energy bounds and average energy for each energy group at the initial 
incident face. In part two, the number of particles in each energy group is 
calculated using the spectral values determined in subroutine FLUXEQ. If 
an integral type spectrum is to be used, the N(>E) is evaluated at the 
boundary energies of each group and the number of particles in each group is 
obtained by differencing the successive values at the boundary energies. The 
differential dN/dE is evaluated at the average energy for each group, and 
the number of particles in each group is computed by multiplying this value 
of dN/dE by the difference in boundary energies for that group. 

The evaporation neutron dose is computed in subroutine EVNEDO. The 
first section contains the input statements for the materials which comprise 
the shield to be analyzed and the computation of the angles to be used for 
the dose calculations . The second part computes and saves two summations 
for each angle and SX sub interval. The first is the summation of the 
product-quotient Hr/P, while the second is the product KSr, with an adjust- 
ment multiplier l sometimes being included to account for the variation of 
materials in the makeup of the slab. Section three computes the evaporation 
neutron source terms for each SX subinterval and also computes the evap- 
oration neutron dose at each print bound based on the above summations. 

Subroutine PROPTY controls the input of data for each material and the 
flux to dose conversion factor tables from the magnetic tape to the computing 
machine. The data tape is mounted on logical tape unit 3 and the mate- 
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rial data is then transferred to disk storage to decrease access time, if 
a disk is available, or to logical unit 4 if a disk is not available. This 
transferring of data reduces the possibility of destroying the master data 
tape in the process of program execution. The flux to dose conversion 
factor tables are transferred directly into computer storage. The second 
section of PROPTY controls the transfer of material data into computer 
storage and initializes the various subroutines involved for the necessary 
ccarstanticand' type oof . int.^p9l^.tion i jtdi3blapp.liddV't'ff' : ER,erdat,a.;;' If.- a call::., 
is made Tor a material which is not available the program will stop and 
print out an appropriate error message. 

The remaining subroutines, XS, LAGRNG, YIELDS, RANGE, CA9NRG, DOSEK, 
and SORT are all straight forward and require no further explanation. 


INPUT DATA FOR LPSC 

Most of the data used by LPSC is on a data tape. This includes the 
range energy tables, secondary particle yield functions, energy of second- 
ary particles, cross sections, mass stopping power, flux to dose conversion 
factors and an RBE table. 

The following data is input on cards: 

NOCDS - the number of identification or comments cards printed at the begin- 
ning of each set of output data. A minimum of one (l) card is required, 
even if blank, and "a maximum of 09 is permissible. Format ,( l£ ) . 


12 


73 80 

XX 


NOJD 


CARD - The comments or identification card(s) to be printed; a total of 
NOCDS required. If print position one (l) on output is for carriage con- 
trol, include appropriate control character in card column one (l) of com- 
ments card. Format (12A6). 


III 

2 

72 

73 80 

■ 

Comments or identification 


ID 


SPBND - The minimum energy in MeV of incident protons which produce second- 
ary particles. 

SNBND - The minimum energy in MeV of incident neutrons which produce sec- 
ondary particles. 
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PDSBND - The minimum energy in MeV for computing proton dose due to ioniza- 
tion. 

NDSBND -r.The minimum energy in MeV for computing cascade neutron dose. 

BNDLOW - The minimum energy in MeV of the initial incident primary proton 
spectrum. 

KNTRP - A control governing the various generations of protons to be cal- 
culated. If KNTRP * 1, primary protons only will be calculated; * 2, 
primary and first generation secondary protons and evaporation neutrons 
produced by primary protons will be calculated; “3, primary and all gen- 
erations of secondary protons and evaporation neutrons produced by all 
protons . 

KNTRN - The control for the generations of neutrons to be calculated. For 
KNTRN = 1, first generation cascade neutrons are calculated; “2, all 
generations of cascade neutrons and evaporation neutrons from cascade 
neutrons are calculated. 

SOFENO - Material number for dE/dX (of receiver) table to be used in dose 
calculations . 


SOFENO 

Material 

1 

2 

Water 

Tissue 


For SOFENO equal to any other number, the program will stop and print 
an appropriate error message. Format (5F6. 0,314). 


1 6 

7 12 

13 18 

19 24 

25 30 

31 

34 

35 38 

39 42 


73 80 

XX.X 

XX.X 

XX.X 

XX.X 

XX.X 


X 

X 

X 


Limit 1 


KOSW - Branching controls for various calculations throughout the program. 
Provision has been made for 36 such controls but not all of them are used. 

The following list describes the effect when the control is set equal to 1 or 
2; for any other number some type of error is likely to occur. (The number 
in parentheses refers to a specific card column.) 

KOSW (l) =« 1, omit table of neutron source terms; 

=2, print table of neutron source terms; 

KOSW (3) m 1, omit table of initial incident energy group bounds, delta energy 
and average energy for each group, value of dN/dE at the average 
energy or N(>E) at energy group bounds, and the number of protons 
in each group; 

=2, print the above table of initial values; 
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KOSW (5) = 1, calculate the dose at each shield print hound; 

= 2, calculate the dose of the initial incident primary proton 
spectrum only without performing any further shield calcula- 
tions ; 

KOSW (7) * 1, energy group bounds for cascade neutron spectrum is the same 
as that for initial proton energy group bounds; 

® 2, energy group bounds for cascade neutron distribution are input 
data; 

KOSW (9) = 1, omit tables of energies of proton particles incident at all 
6X layers and associated proton cross-sections; 

-2, print tables of energies of protons and associated cross- 
sections. (This data can be used for checking purposes.); 

KOSW (ll) = 1, omit tables of primary and secondary proton and cascade 

neutron spectrum source terms at intermediate print bounds 
and exit face; 

=2, print tables of all spectrum source terms at intermediate 
print bounds and exit face; 

KOSW (13) = 1, initial incident proton spectra data is not time integrated; 

= 2, calculations are based on time integrated initial incident 
proton spectra data; 

KOSW (15) =* 1, construct additional energy groups at intermediate print 
bounds using equal increments; 

=2, construct additional energy groups at intermediate print 
bounds using variably spaced increments; 


KOSW (17) s* 1, omit table of total proton and cascade neutron flux terms; 

-2, print table of total proton and cascade neutron flux terms; 

KOSW (19) « 1, initial incident energy group boundaries are calculated from 
the energy group bounds at the exit face and intermediate 
print bounds; 

« 2, initial incident energy group boundaries are input data. 
(This feature is useful for a monoenergetic case) ; 


KOSW (21) = 1, calculation of proton dose factors due to nuclear interaction 
is bypassed; 

=2, proton dose factors due to nuclear interaction are calculated. 


This version of the program is equipped to handle this type of dose computa- 
tion; however, at present there is no reliable data for these calculations 
and KOSW (21) should be set equal to 1. If reliable data is forthcoming it 
may be used in place of the table of zeros now in the program. 
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For the following control switches (card columns 25 through 36), the 
program uses linear interpolation on the raw data when the branching control 
equals 1; and for the control switch equal to 2 the program replaces the raw 
data by the logarithms (base 10) of the data and then uses linear interpo- 
lation on the logarithms. The interpolation with logarithms option was in- 
cluded to increase the accuracy of the calculations, however the use of this 
option will also increase the computing time necessary for any problem. 


KOSW (25), 
KOSW (26), 
KOSW (27), 
KOSW (28), 
KOSW (29), 
KOSW (30), 
KOSW (31), 
KOSW (32), 
KOSW (33), 
KOSW (34), 
KOSW (35), 
KOSW (36), 


range -energy; 

proton cross-section; 

neutron cross-section; 

emitted yields for protons bombarding; 

emitted yields for neutrons bombarding; 

energy of cascade particles for protons bombarding; 

energy of cascade particles for neutrons bombarding; 

dE/dX mass stopping power; 

relative biological effectiveness (RBE); 

cascade neutron flux to dose factors for rad or rad/hr units; 

cascade neutron flux to dose factors for rem or rem/hr units; 

proton (nuclear interaction) flux to dose factors. (This table 
presently contains zeros.) Set KOSW (36) equal to 1. 


FORMAT (3611) 


1 2+ 36 


73 80 

xxxxxxxxxxi XXXXXXXXXXXl 


KOSW 


EQX - The number of shield thicknesses or print bounds at which data is 
desired. The maximum number permissable is 20. 


NOAECt - The number of angles used in computing the evaporation neutron 
dose. This must be an integer in the range from 1 to 10. 


Format (214) 


14 5 8 

XX XX 


73 80 

Limit 2 
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The sevan pieces of data which follow are needed for each shield thickness 
at which data is to be printed - a total of NOX cards. 

X( J) - The shield thickness in g/cm^ at the print bound; 

NQD2X ( J ) - The number of SX increments between the and (J-l) th print 
bounds. If J = 1, the (j-l)^* 1 print bound is assumed to be at X = 0. 

The total number of 5X increments for the entire shield must be < 200. 

FROPNO(j) - The property number for the material between the 1 th and 
( J-i)tb pyint bounds. This version contains data for the following mate- 
rials : 


Material 

PROPNO 

Hydrogen 

2 

Beryllium 

104 

Carbon 

102 

Nitrogen 

106 

Oxygen 

107 

Aluminum 

101 

Titanium 

108 

Iron 

103 

Tungsten 

109 

Lead 

105 

Uranium 

110 

Water 

1 

Polyethylene 

3 


If a PR0PN0( J) is set equal to any other number, the program will stop and 
print an appropriate error message. 

P( J) - The density of the material in g/cm^ 

H( J) - The ratio of the hydrogen density in the material to that in 
water; 

S( J) - The removal cross-section for the J^* 1 material in cm^/g 

K( J) - The value of K ^ used in the evaporation neutron dose calculations. 
For 2 < Z < 6, K = 1.0; and for Z > 6, k = 0.5. K must be a non-zero 
quantity even though no value for K-^ may be required. 

Format (F6.0, 214, 3F8.0, F6.0) 


1 6 

7 10 

11 14 

15 22 

23 30 

31 38 

39 44 

■ 

73 80 

XXX. XX 

XX 

XXX 

x.xx*x 

x.xxxx 

x.xxxxx 

x.x 

■ 
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If K0SW(l9) = 1, the following data is required to establish the energy- 
group bounds at the exit face and intermediate print bounds. 

1.) NOINTS(l) - The number of groups of equal energy increments at the 
exit face, a maximum of 25 groups permissible. 


Format (13) 


1 3 
xx 


73 80 

1 or 
NOINT 2 


Then NOINTS(l) pairs of numbers defined by: 

2. ) E0MAX( J ,l) - The maximum energy in MeV of the group of equal 
energy increments at the exit face. For example, if the fourth group of 
intervals has a maximum energy of 100 MeV and the fifth group 200 MeV, then 
E0MAX(4,l) = 100 and E0MAX(5,l) « 200. 

3. ) N0HJCR(J,1) - The number of energy increments in the Jfch group at 
the exit face, i.e., if the fifth group has 5 increments, then R0INCR(5,l) «* 5 
and the increment size 5E = (200-100) /5 ** 20 MeV. The total number of energy 
increments at the exit face or any intermediate print bound cannot exceed 300. 
Use as many cards as required with 6 pairs of number per card. 


Format (F8.0, 14, F8.0, 14, F8.0, 14, F8.0, 14, F8.0, If, F8,0, 14) 



9 12 

13 20 

21 24 

25 32 

33 36 

37 44 

45 48 

49 56 

57 59 

61 68 

.ms 

XX 

xxxx.x 

XX 

xxxx.x 

xx 

xxxx.x 

XX 

xxxx.x 

XX 

xxxx.x 


CD 

CD 

73 

80 


1 

or 

XX 

EBJT2 



For K0SW ( 15 ) = 2, a set of data similar to l) through 3) is needed which 
applies to the initial face and all intermediate print out boundaries i.e., 

4) N0IETS(2), 5) E0MAX(J,2), and 6) N0INCR ( J, 2 ) . Whereas at the exit face 
the energy grouping applies to the whole energy grid, at the initial face and 
print bounds these groups are used to determine the energy grid between 
E = E 0 and E = 0 only, where E 0 is the energy on the initial face or 

print bound which degrades to E * 0 at the (K + l)^* 1 intermediate print 
bound or exit face. The remainder of the energy grid at the initial face or 

print bound is calculated from the energy grid at the (K + l)^* 1 print 

bound or exit face. If K0SW(l5) = 1, data 4) through 6) are replaced by 

4'. ) NDE - the number of equal energy increments into which the range E = 0 















24 


to E = E 0 will be divided at the initial face and intermediate print out 
bounds, i.e., in the range (0,E o ) there will be NDE increments of size 
5E = Eq/NDE. Format (13). 


1 3 


73 80 


xx 


NDE 


7.) E3MAX - The maximum energy bound in MeV at the initial incident face. 
Format (E12.5) 



The monoenergetic case or any case for which the energy grid at the 
initial incident face is to be input can be computed by setting K0SW(l9) = 2, 
omitting data words 1.) through 7.) above, and inserting the following 
information: 

.'.'A KBl - The : htunber of energy bounds at the initial incident face; 


Format (13) 


1 3 


73 80 

XX 


KEI 


El - The energy bounds for each group, a total of KEI. The maximum number 
permitted is 300. Format (8F9.0) 


1 9 

10 18 

19 27 


37 45 

46 54 

55 63 

64 72 

mm 

xxxx.x 

XXXXuX 

xxxx.x 


xxxx.x 

xxxx.x 

xxxx.x 

xxxx.x 

1 


The spectrum simulation is controlled by two or more cards as follows: 

MOVE - A control for the type of spectrum to be used. If r(E) =* dN/dE, 
MOVE « Z, and for T(E) - N(>E), MOVE » 1. 

EQNO - The number of the equation or method describing the proton spectrum 
at the initial incident face. Seven different forms, numbered 1 through 7, 
are available and described below. 

TITLE - An identification for the title or name of the spectrum. A maximum 
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of 66 characters, including blank spaces, are permitted for this identification. 
Format (213, 11^6) 


1 3 

4 6 

7 

72 

73 80 

X 

X 

Name of spectrum 


CONTRL 1 


If EQNO - 1, T(E) » AE" B 

The coefficients A and B in the equation are input. Format (2E12.5) 


1 12 
±x.xxxxxE±yy 


13 24 

±x.xxxxxE±yy 


73 80 

CONTRL 2 


If EQNO - 2, F(P) - N(>P) -A exp(-P/P Q ) 

MOVE = 1 on CONTROL 1 and the coefficients A and P Q replace A and B 
on CONTROL 2. Format (2E12.S) ™ 

IF EQNO * 3, a table of values is to be read as input. 


In place of card CONTRL 2, read the following: 

NOENTS - The number of entries in the table on card FLUX 3,1. A maximum of 
100 entries in the table is permitted. Format (14) 


1 4 


73 80 

XX 


FLUX 3,1 


This is followed by the table in Format (F8.0, E10.3, F8.0, E10.3, F8.0, E10.3, 
F8.0, E10.3). These represent pairs of numbers, the energy' - FF.ETF! and the 
corresponding spectra value FROTS at this energy. 


1 8 

9 18 

19 26 

27 36 

37 44 

45 54 

55 62 

XXX. XX 

+x.xxxE±yy 

XXX. XX 

+x.xxxE±yy 

XXX. XX 

+x.xxxE±yy 

XXX. XX 


63 72 


+x . xxxEtyy 

illBSi 


Continue with as many cards as required. 
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4 

■2 


-E 1 ' 1 and B(E) = 


-(.=)■■= 4 

where A(E) 

i«l 

and b are input. 

a^, i = 1, 2, 3, 4. Format (4E12.5) 


If EQNO = 4, T(E) = A(E) exp(-B(E)) 

4 


Z ^ El "- 


The coefficients a 


i=l 


1 12 
±x . xxxxxElyy 



E 


73 $0 

FLUX 4A 

bjy 4 = 1* 2, 

3, 4. Format 

(4E12.5) 

1 12 
±x.xxxxxE±yy 

13 24 

±x . xxxxoflBkyy 

25 36 

tx.xxxxxEiyy 

37 48 

±x . xxxxxE+yy 


73 80 

FLUX 4B 


If EQNO = 5, log 10 r(E) = ^ a ± : 

The a-j^ are input as in card FLUX 4A. 


E 


1-1 


i-1 


If EQNO = 6, log 10 r(E) 


a i (log 1Q E) 


i-1 


i=l 


The a^ are input as in card FLUX 4A. 

If EQNO = 7, r(P) - || - A e' P / P ° 

MOVE = 2 on CONTRL 1 and the input A. and P Q are in FORMAT (2E12.5) as on 
CORTRL 2. 

For. EQNO = 1, 3, 4, 5, 6, MOVE may be either 1 or 2 depending on the 
data to be used. 


The last set of cards construct the cascade neutron energy table. The 
first card in this set gives the number of entries, NONUBD, in the table. 
Format (13). 
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3 

xx 


(' 

73 80 


N0NUBD 


This is followed by the bounds, NUENBD, of the energy groups in descending 
order ... 100 MeV, 80 MeV, 60 MeV, ... 0 MeV. Format (7F10.0). 


1 10 

11 20 

21 30 

31 60 

73 80 

XXX. X 

XXX. X 

XXX. X 

• • • 

NUENBD 


Use as many cards as required. 


TAPFIX Program 

Program TAPFIX is a FORTRAN IV computer program which generates from 
tables of data punched in cards a data tape to be used with the Lewis Proton 
Shielding Code (LPSC). The tape to be generated is mounted on logical unit 3. 
The first half of TAPFIX deals with the property data for various materials 
and the second part controls the tables of flux to dose conversion factors. 
Following is a list, in sequential order, of the data and formats necessary 
to generate a tape which is compatible with LPSC. 

NOMAT 1 - the number of various materials for which shielding data is 
available . 

NOMAT 2 - the number of receiver materials for which tables of . dE/dx (to be 
used in dose calculations) are available in this program. 


FORMAT (214) 


1 4 

5 

8 


73 80 

XX 


XX 


LIMITS 


Following this card should be NOMAT 1 sets of data, one for each material. 

MATNO - the number assigned to each material (see list on page 22); hydro- 
genous materials are in the range 1 through 49, nonhydrogenous materials are 
assigned numbers greater than 100. 

N0FC0M - the number of elements comprising a given material. All the materials 
in LPSC are single element materials except water and polyethylene which 
contain two elements. The program is coded to allow compounds with up to 
four elements. The table entries for secondary yields of compounds are dif- 
ferent from the yields for elements in the following way. For a compound. 
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the yields for each element are multiplied by the ratio of the inelastic 
cross section of the i^* 1 element to the total inelastic cross section for 
the compound, thus the fraction of incident particles that collide with the 
i^* 1 nucleous produce yields from the i^h nucleous etc. This requires sec- 
ondary yields and secondary energies for each’ element of bach compound-iV - 
If the compound is hydrogenous, the hydrogen data must be placed first in 
each group of tables. The cross section tables for neutrons and protons 
bombarding are for the compound or element which ever applies. 

GMWT - gram molecular weight of the particular material. 

LI - length of the range-energy table. 

L2 - length of the energy of cascade particle table - protons bombarding. 

L3 - length of the energy of cascade particle table - neutrons bombarding. 

L4 - length of the emitted yield table - protons bombarding. 

L5 - length of the emitted yield table - neutrons bombarding. 

L6 - length of the proton cross-section table. 

L7 - length of the neutron cross-section table. 

FORMAT (214, E13.6, 714) 


1 4 

5 8 

9 21 

22 25 



mm 

gum 


46 49 

XXX 

X 

±x.xxxxxxE±YY 

XX 


XX 



XX 

XX 


73 

LKDj 


ENERGY - energy grid in MeV for range -energy table. 

RANGE - associated range for each entry in ENERGY grid. (LI pairs of numbers) 
FORMAT (8E9.3) 


1 9 

10 18 

19 27 

28 36 

37 45 


64 72 

73 ' 8< 9 

±x.xxxE±x 

±x.xxxE±x 

±x.xxxE±x 

...... ml — 

±x.xxxE±x 

±x.xxxE±x 

• # •_ 

±x.xxxE±x 

±x.xxxE±x 


ENRGRP - energy in MeV of bombarding proton for cascade particle energy table. 
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EPRPR y associated energy in MeV of cascade proton produced. 

EFRMJ - associated energy in MeV of cascade neutron produced (L2 triads of 
numbers) . 

FORMAT (9E8.2) 


1 8 

9 16 

17 24 

25 32 

33 40 

41 48 

49 756 

57 ' 64 

±x.xxE±Y 

±x.*xE±Y 

±x.xxE±Y 

±x.xxE±Y 

±x.xxE±Y 

±x.xxE+Y 


±x.xxE±Y 


65 72 

73 

80 

ix.xxEiY 




EERGEU - energy in MeV of bombarding neutron for cascade particle energy 
table . 

ENUPR - associated energy in MeV of cascade proton produced. 

EHUNU - associated energy in MeV of cascade neutron produced. (L3 triads 
of numbers) 

FORMAT (Same as for energy of cascade particle table - protons bombarding) 

EKERPR - energy in MeV of bombarding proton for emitted yield table; 

YPRCP - associated cascade proton yield per inelastic event; 

YPRCN - associated cascade neutron yield per inelastic event; 

YPREM - associated evaporation neutron yield per inelastic event; (L4 groups 
of numbers) 

FORMAT (8F9.0) 


1 9 

H 

o 

19 27 

28 36 

LO 

r*- 

to 

46 54 

55 63 

64 72 

73 

80 

XXX. X 

x.xxx 

X.XXX 

x.xxx 

XXX. X 

x.xxx 

x.xxx 

x.xxx 




EEERNU - energy in MeV of > bombarding neutron for emitted -yield table; 

YMJCP - associated cascade proton yield per inelastic event; 

YNUCN - associated cascade neutron yield per inelastic event; 

YMJEN - associated evaporation neutron yield per inelastic event; (L5 groups 
of numbers) 
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FORMAT (Same as for yield tables with protons bombarding) 

EBOMP - energy in MeV of proton particle for cross-section table; 

XSPR - associated proton cross section in millibarns; (L6 pairs of 
numbers) 

FORMAT (1QF7.0) 


1 7 

8 14 

15 21 

22 28 

2977 35 

36 42 

43 49 

50 56 

57 63 

64 70 

XXX. X 

XXX. X 

xxx.x 

xxx.x 

xxx.x 

xxx.x 

xxx.x 

xxx.x 

xxx.x 

xxx.x 


73 80 


EBOMN - energy in MeV of cascade neutron particle for neutron cross section 
table . 

XSMJ - associated cascade neutron cross section in millibarns (L7 pairs of 
numbers). 

FORMAT (Same as for proton cross section table) 

This concludes the liBt of property data for each material. The follow- 
ing data are flux to dose conversion factors for protons and neutrons. 

L8 - length of the relative biological effectiveness table (RBE) 

L9_ - length of neutron flux to dose factors table for rad or rad/hr units 

L10 - length of neutron flux to dose factors table for rem or rem/hr units 

Lll - length of proton (nuclear interaction) flux to dose factors table 

FORMAT (414) 


1 4 

5 8 

9 12 

13 16 


73 80 

XX 

XX 

XX 

XX 




RBEERG - energy value in MeV for the RBE table 

RBE - associated relative biological effectiveness (L8 pairs of numbers) 
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FORMAT (F7.0, E9.3, F7.0, E9.3, . . . ) 


PH 


iim^i 

24 32 


40 48 

49 55 

56 64 





±x.xxxE±Y 

mi 

±x.xxxE±Y 

XXX. X 

±x . xxxE±Y 



73 80 


K1NRG - energy value in MeV for neutron flux to dose conversion table in 
rad or rad/hr units. 

K1 - associated neutron flux to dose conversion factor in (rad/hr)/ 

/neutrons / cm 2 sec) (L9 pairs of numbers). 

FORMAT (Same as for RBE table). 

K2NRG - energy value in MeV for neutron flux to dose conversion table in 
rem or rem/hr units. 

K2 - associated neutron flux to dose conversion factortin 
Xneutron/cm 2 sec) (L10 pairs of numbers). 

FORMAT (Same as for RBE table) 

KNRG - energy value in MeV for proton (nuclear interaction) flux to dose 
conversion table. 

KK - associated proton flux to dose conversion factor ( Lll pairs of numbers) . 

FORMAT (Same as for RBE table). This last table has zero value entries in 
this version of the code due to the poor data available. 

The last data entered on the tape is NOMAT 2 sets of dE/dx tables; 
one for each receiver material chosen. 



MATNO 2 - the receiver material number 
(See table on page 19) 

L12 - length of the dE/dx table. 
FORMAT (214) 

1 4 1 5 8 ~ 


73 


80 


xx 


XX 
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EXEHRG - energy value in MeV for dE/dx dose conversion table. 

DEDX - associated dE/dx conversion value in MeV cm^/g (L12 pairs of num- 
bers) 

FORMAT (Same as for RBE table) 


Sample Problem 


A sample problem constructed for instructional purposes follows. The 
problem was to calculate spectra and doses at three print boundaries for a water 
shield 30 g/cm 2 thick. The print bounds were chosen at 10 g/cm 2 , 20 g/cm^ and 
30 g/cm^. A time integrated type proton spectrum was used where 


N(>E) = 7.45X10 12 E“ 2,12 protons/cm 2 

when the values at 10 g/cm 2 are calculated the program assumes that the shield 
consists of 10 g/cm^ only and similarly for the calculations at 20 and 30 g/cm^. 


The primary protons and all generations of cascade protons, cascade neu- 
trons, and evaporation neutron are calculated. The isotropic emission of 
evaporation neutrons is simulated using five angles. 


Table I contains the input data for the sample problem. The first card 
indicates the number of comments cards which follow. The card labeled LIMIT 1 
contains the cut-off energies for the production of secondary particles by 
protons, secondary particles by neutrons, proton dose calculations, neutron 
dose calculations, incident proton energy, and two control numbers used to 
control the calculation of secondary particles and a receiver material number 
respectively. The next card labeled KOS'W controls the print of data and the 
manner of interpolation in the various data tables. The next card containing 
a 3 and 5 indicate the number of print bounds and the number of angles used 
in the evaporation neutron calculation, respectively. The number of angles 
used in the evaporation neutron calculation is < 10. Since the running time 
of the program is angle dependent, five angles are recommended. There are 
cases where the evaporation neutron dose is less than 1 or 2 percent of the 
total dose. For these cases one may want to use one angle which will result 
in the evaporation neutrons running faster. The next 3 cards contain the 
print bounds, the number of 5X contained in AX, the material' property num- 
ber, the material density, the hydrogen ratio, the removal cross section 
(oxygen only) and the value of K. The next card containing the number 6 
indicates the number of AE boundaries where 6E may change for the exit 
face. These represent nonzero energies at the exit face. The next card 
indicates that there exists 6 energies between 0 and 6 MeV, 2 energies be- 
tween 6 and 10 MeV, 5 energies between 10 and 60 MeV, etc. The next set 
of data constructs the exit energy bounds at intermediate print boundaries . 

The next card labeled E Max is the maximum energy of incident protons to 
be considered. The next two cards labeled FLUX 1 define the proton spec- 


E-3080 
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TABLE I. - INPUT DATA FOR SAMPLE PROBLEM 


2 

SAMPLE PROBLEM 
WATER SHIELD 


10. 10.0 2.0 

0.0 20.00 

3 

2 2 



LIMIT 1 

221212722 1 

1 222222 

2 2 2 2 0 1 




KOSW 

3 5 






LIMIT 2 

10. 20 11.0 

1.0 

.033 

1 .0 



MAT 

20. 10 1 1.0 

1.0 

.033 

1 .0 




30. 10 1 1.0 

1.0 

.033 

1 .0 




6 






NOINT 1 

6. 6 10. 

2 6 0; 

5 

100. 

2 60 C. 

5 1000. 

2 EINT 1 

7 






NOINT 2 

6 . 6 1 0 . 

2 20. 

4 

40. 

10 100. 

20 6C0. 

3 EINT 2 

100C. 2 






EINT 2 

+1 . 000CCE+P3 






EIMAX 

1 1N=A*E**(-B) 






CONTRL 1 

+ 7.4547CE + 12 J -2. 1200 

OE+P 





CONTRL 2 

17 






NONUBD 

1000. 400 . 

3 00. 

■>ro. 

ICC. 

3 0. 

60. 

NUENBD 

50. 40. 

30. 

20. 

I '«/ • 

3 . 

6 . 

NUENBD 

4 . 2. 

0. 





NUENBD 


TABLE II. - INPUT DATA PRINTED WITH OUTPUT 


SAMPLE PROBLEM 
WATER SHIELD 


X »M IN 

X.MAX 

NUMBER OF 

MAT E R I AL 


density 

HYDROGEN 

REMOVAL XSFCT 

( GM /CM** 2 ) 

INCREMENTS 

NUMBFR 

< 

GM/CM**3 ) 

RAT I 0 

( CM **2 /GM ) 

0 . 

10. CO 

20 

l 


1.0000 1 

1.000 

0.03 30 

10. CO 

20.00 

10 

1 


1.00000 

1.000 

0.0330 

20.00 

30.00 

10 

l 


1.00000 

1.000 

0 . 0 3 30 

NUMBER 

CF ANGLES IN EVAPORATION 

NEUTRON DOSE 

CALCJL4TI ONS 

= 5 


ENERGY 

CUT-OFF 

LEVEL FOR INITIAL 

INC I DENT 

PROTON FLUX = 

20.00 MFV. 


ENERGY 

CUT-OFF 

LEVEL FOR CALCULATING SECONDARY PROTONS = 

10.00 M F V . 


ENERGY 

CUT-OFF 

LEVEL FOR CALCULATING CASCADE 

NEUTRONS = 

10.00 MEV. 


ENERGY 

CUT-OFF 

LEVEL FOR PRGTCN 

DCSE DUE 

TO 

I UN I ZA T I ON = 

2.00 3EV. 



ENERGY CUT-OFF LEVEL FOR NEUTRON DOSE CALCULATIONS = 0. MEV. 

® = A*E**(-B) with A = 7.45470E 12 ANO B = 2.12DOOF 00 


N=A*E** ( -B ) 


TABLE III. - PROTON ENERGY AND SPECTRUM DATA ON INCIDENT FACE OF SHIELD 
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trum. The next card with 17 on it indicates the number of energy boundaries 
used for constructing the cascade neutron energy groups. The following data 
are the cascade neutron energy boundaries. 

Table II is a printout of some of the input data. The first two com- 
ments cards describe the purpose and type of shield considered. The three 
lines of data under the title block show the upper and lower bounds of Ax 
(labeled x,min and x,max, respectively, in g/cm^). Also shown are the 
number of increments in each Ax, the material number (this is a call number 
for getting data from the tape), the density, the hydrogen ratio, and the 
removal cross section. This data was kept on cards so it could be readily 
changed. Below this data is indicated the number of angles used in the 
evaporation neutron calculation. The next five lines indicate the ..cut-off 
energy for the indicated calculations. Below the cut-off energies is the 
spectrum equation and the coefficients used. 

Table III is a print of the incident proton spectrum and the energy 
boundaries and energy groups used at the incident face. These energies 
were calculated from input on the exit face. The break between 22 and 23 
(left column) shows the maximum energy 2.16621E02 which is X) at the exit 
face. Also 2.16612E02 is the first energy group that does not penetrate 
to the exit face. Similar breaks occur for each print bound. Table III 
shows six columns of data. The titles should be self explanatory. The 
column' on the left without a title is an index counter showing how many 
entries there are in the table. Reading from left to right the second 
column shows the incident proton energy bounds at the incident face of the 
shield. The next column is N(X) for the indicated energy bounds in column 
two. Column four is the difference of successive values in column two. 
Column five is the midpoint energy of column two. The last column is the 
difference of successive values of column three and represents protons/cm 2 
or protons /cm? s ec. in group i. 

If a differential type spectrum were used the column headings for 
table III would have been: number of entries, proton energy bounds on the 

incident face in MeV, 5E in MeV, average energy in MeV, the differential 
spectrum in protons / cm^ MeV or protons/cm 2 MeV sec and the (cLN/dE)SE in 
protons/cm^ or protons/cm^ sec in each group. 

Table IV contains the proton spectrum data for the print bound 
x = 30. g/ cm 2 . The data for the print bounds x = 10 and 20. g/cm 2 were 
delfeted from the report. As noted on the incident face data, table III 
22 energy groups penetrate the shield. The remainder of the 139 incident 
energies were stopped by the shield. The zero data between group 24 and 
group 138 was not included in the report. The column on the left indicates 
the number of entries in the table. Reading from left to right, the second 
column indicates the midpoint energy of each proton energy group in MeV at 
the indicated thickness. The entries in the third column are the number of 
primary protons in each of the indicated energy groups. These are in 
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protons/cm? for time integrated spectra or protons/cm^ sec for spectra that 
has not been time integrated as indicated at the top of the table. The 
following two columns contain the cascade secondary protons at thickness x. 
The column labeled FIRST GEN contains the first generation cascade protons 
at x. The next column labeled HIGHER GEN contains the second and higher 
generation cascade protons at x. The last two columns show the same param- 
eters as the previous two columns, however, these are from the last small 
Sx only before reaching x. At the bottom of this table is indicated the 
total number of particles in the indicated columns and the rad and rem dose 
respectively from these components. 

Table V is read similarly to table IV only these data are for cascade 
neutrons at x and from the last small 6x only. These data are also for 
first generation and second and higher generation particles as indicated. 

Table VI contains the dose in rad per flare for the various radiation 
components calculated. The columns are labeled to be self explanatory. 

Table VII is read the same as table VI only the units are rem per 
flare. If the spectrum used was in particles/cm^ sec the tables VI and 
VII would be in rad/hr and rem/hr, respectively. 

Table VIII contains the total integrated flux for the various particles 
at the indicated print bounds. The flux of evaporation neutrons was not 
calculated in this table. 

Table IX contains the evaporation and cascade neutron source terms 
for the various layers. The neutron source terms are in neutrons/g or 
neutrons/g sec. The location at these source terms is the midpoint of 
all Sx layers . 


Sample of Nuclear Interaction Data 

The main source of nuclear interaction data used in LPSC which include 
inelastic cross sections, secondary yields, and energy of secondaries (for 
nuclei of mass numbers >1.2) was supplied by H. Bertini of ORNL (see ref. 4). 
This data was supplied for 10 elements ranging from carbon (12) to Uranium 
(238) . The program developed by Bertini for estimating the secondary yields 
and the energies of the secondaries does not give good statistics for nuclei 
containing few nucleons. Carbon was the lowest mass number believed to give 
reasonable statistics. Therefore Bertini* s data was used for mass numbers 
> carbon in this program. The data when plotted at constant energy as a 
function of mass number falls on smooth curves; thus interpolation relative 
to mass number may be accomplished for elements not calculated by Bertini. 

The data supplied in reference 4 is voluminous and it was necessary to 
select from it just the data required for LPSC. A sample of the data re- 
quired for carbon appears in table X. The top set of data is for protons 
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bombarding. Column one is the energy in MeV of the incident particle. Column 
two is the average yield of cascade protons per collision. Column three is 
the average energy in MeV of the emitted protons. Column four is the average 
number of cascade neutrons emitted per collision. Column five is the average 
energy in MeV of the emitted cascade neutrons. Column six is the average 
yield of evaporation neutrons per collision. Column seven is the inelastic 
cross section in millibarns. The bottom set of data is for incident neutrons 
bombarding'. 

The data was all rounded to four significant figures. No claim is made 
that more than one or two are reliable. The data was plotted as a function of 
incident particle energy. The values read from smooth curves drawn through 
the data in table X were used on the data tape. 


Nuclear Interaction Data for Hydrogen 

When high energy nucleons collide with hydrogen nuclei this program (LPSC) 
assumes that scattering of the incident particle occurs. It is further as- 
sumed that the hydrogen nucleus is added to the beam as an additional proton. 

The data in table XI contains the secondary yields, the energies of 
secondaries, and the cross sections for neutrons and protons bombarding 
respectively. This table is similar to the previous table X. 


Source of Data 

Table XII on page 45 was included to show the source of data for proton 
and neutron cross sections, proton range energy, and proton mass stopping 
power. 
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TABLE X. - CARBON 


Yield and energy for protons bombarding 


Incident 
energy 
in MeV 

Average 
yield of 
emitted 
cascade 
protons 
per 

collision 

Average 
energy of 
emitted 
protons 
in 
MeV 

Average 
yield of 
emitted 
cascade 
neutrons 
per 

collision 

Average 
energy of 
emitted 
neutron 
in 
MeV 

Average 
yield of 
evaporation 
neutrons 
per 

collision 

Inelastic 

cross 

section 

in 

millibarns 

25 

0.5806 

9.990 

0.4153 

9.382 

0.02697 

447.9 

50 

*'•.8788 

18.86 

.6338 

17.57 

.1713 

348.9 

100 

1.177 

38.09 

.8039 

33.02 

.2788 

271.6 

150 

1.300 

57.97 

.8669 

52.77 

.2820 

245.5 

200 

1.435 

77.41 

.8507 

69.38 

.3307 

232.3 

250 

1.499 

97.23 

.8900 

82.57 

.3598 

222.5 

300 

1.580 

112.9 

.8864 

103.0 

.3249 

215.4 

350 

1.593 

132.8 

.9325 

114.9 

.3509 

218.0 

400 

1.642 

154.9 

.9250 

124.4 

.3497 

233.4 

Yield and energy for neutrons bombarding 

25 

0.4257 

9.095 

0.5698 

9.865 

0.4466 

443.8 

50 

.6290 

17.16 

.8616 

19.31 

.5125 

353.6 

100 

.7898 

34.62 

1.159 

38.52 

.5848 

266.7 

150 

.8450 

51.97 

1.339 

57.26 

.5526 

236.8 

200 

.9075 

66 . 85 

1.424 

76.58 

.5756 

224.4 

250 

.9157 

80.08 

1.493 

97.46 

.5800 

214.9 

300 

.8780 

102.3 

1.563 

115.2 

.5122 

215.9 

350 

.9038 

116.0 

1.556 

138.2 

.4932 

215.7 

400 

.9170 

125.1 

1.677 

151.1 

1 

.5301 

223.8 
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TABLE XII. - REFERENCE TABLE FOR CROSS SECTIONS, RANGE ENERGY, AND MASS 


stopping; power data 


Material 

Mass stopping power 
and range energy 

Low energy neutron 
cross sections 

High energy proton and 
neutron cross sections 

Hydrogen 

Private communica- 

Reference 9 

References 9, 10 


tion, C 

!.W. Hill 




Lockheed, Georgia 



Beryllium 



No secondaries 

Proton attenuation by 




being generated 

ionization only 

Carbon 



Reference 11, 15 

Reference 4 

Nitrogen 



Reference 12 

Interpolated data from 





reference 4 

Oxygen 



Reference 12, 15 

Reference 4 

Aluminum 



Reference 12 

Reference 4 

Titanium* 



Reference 12 

Interpolated data from 





reference 4 

Iron 



Reference 13, 15 

Interpolated data from 





reference 4 

Tungsten 



Reference 14, 16 

Reference 4 

Lead 



Reference 14, 16 

Reference 4 

Uranium 



Reference 15 

Reference 4 

Water 



Reference 9, 12 

References 4, 9, 10 

Polyethylene 

> 

f 

References 9, 10, 

References 4, 9, 10 
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The Lockheed range energy data was interpolated for titanium. 
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APPENDIX A 


o 

oo 

0 

to 

1 


w 


LPSC PROGRAM LISTING 

IBFTC P I SR DECK, LIST 

THE LEWIS PROTON SHIELDING CODE ( LPSC) 

PROTON INDUCED SECONDARY RADIATION. 

THIS PROGRAM CALCULATES THE PRIMARY PROTON, CASCACE PROTON, 

CASCADF NEUTRON, AND EVAPORATION NEUTRON DOSES IN RAD AND REM 


FOR A BEAM OF PROTONS 

INCIDENT NORMAL 

TO A SLAB 

COMMON 

D2XI20) 

♦ 

MAX 

, EO ( 300, 20 ) 

COMMON 

El (300) 

♦ 

DEI (300) 

, EIBAR ( 300 ) 

COMMON 

OP ( 300 ) 


OPPRM (300 ) 

, N0D2XI20) 

COMMON 

XI 20) 

* 

NOX 

, ENTOTSf 200 ) 

COMMON 

DXI20) 

f 

PROPNO! 20 ) 

. Cl 

COMMON 

PDSBND 

* 

NDSBND 

, BNDLOW 

COMMON 



ENRG( 100) 

, RNG(IOO) 

COMMON 

EBOMBP (25.4) 

* 

EBOMBN( 25*4) 

, CPSP (25,4) 

COMMON 

CPSN I 25 , A ) 

♦ 

CNSPC 25,4) 

, CNSN ( 25,4) 

COMMON 

ENSP (25,4) 

♦ 

ENSN( 25,4) 

♦ EBOMP (25,4) 

COMMON 

EBOMN ( 25 » A ) 

* 

EPROP (25 » A ) 

, EPRON (25,4) 

COMMON 

ENEUP (25*4) 

9 

ENEUN (25,4) 

, ENRGP ( 25 ) 

COMMON 

FNRGN ( 75 ) 

9 

XSMBP (25) 

, XSMBN(75) 

COMMON 

SNRG( 100) 

9 

RBFNRGt 20 ) 

, C1NRG (40 ) 

COMMON 

C2NRGI40) 

9 

CNRG( 2) 

, SOFF(IOO) 

COMMON 

RBF ( 20 ) 

9 

C0NKK40) 

, C0NK2 ( 40 ) 

COMMON 

CONK ( 2 ) 

9 

LENGTH! B ) 

♦ GMWT 

COMMON 

LSOFE 

9 

LRBE 

, LK1 

COMMON 

LK2 

9 

LK 

, NOFCOM 

COMMON 

MOVE 

9 

KOSW (36) 

, KONSWT 

EQUIVALENCE (KOSW! 1) 

.KOSWl ) »{KOSW( 5 ) .K0SW5 ),(KOSW( 7),K0SW7 


A ( KOSW ( 9),K0SW9 ) , ( KOSWI 1 1 ) ,K0SW1 1 ) , ( KOSWI 17 ) ,K0SW13 ) , 

P ( KOSW (17) ,K0SW17) , ( KOSWI 71 ) ,K0SW21 I 

DIMENSION CARD(12),NUFNBD(300),NUENRG(300) ,S I GCNP ( 300 , 2 ) , 

APHI JIf 300,3) , EIBNDS { 300 ) ,AVGNRG(300) »NEUTXS(300) »DIST{300) « 
BSIGNEL(30C) ,UCNEUT ( 300 , 2 I , 

CUCPPIM ( 300 ) ,UCSEC (300,2 ) ,CPPf300,2J ,CNP(300,2) ,PROTS ( 300 * 5 ) 

D » NEWTS ( 300,4) , TOTALS ( 20,9) ,PDOSF( 20 , 5 , 3 ) .PdOsRM ( 20 , 5 ,3 ) ,ND0SE(20,4 
E) »NDOSRM( 20,4 ) , FVNUDS (20) ,EVNDRM( 20) , ANSI 3) , CAS ANSI 2) 

DIMENSION TSPI3) *TP ( 3 ) , TD ( 3 ) , XM ID ( 200 ) .CNTOTS I 2 00 * 2 ) ,SOTEEN I 200 ) 


A , MI DNRG I 300 ) 

REAL NUENBD*NUENRG,NEUTXS, NEWTS, NDOSE , NDOSRM, KOEE 

REAL NDSBND,IPDRAD»IPDREM, IPFLTO 

REAL MI DNRG 

INTFGFR SOFENO 

INTEGER PROPNO 

EQUIVALENCE (PROTStl ) ,PHIJI 1 1 ) ) , I PROTS 1 901 ) ,CPP(1) ) , I NEWTS 1 1 ) , 
ASIGCNP 1 1 1 ) . (NFWTS(601 ) ,CNP( 1 ) ) 

I ANS 1 1 ) »CPS ) * IANS 1 2 ) 

I ENRGYP ,CASANS(l)!*ll 
I EO I 1 ) ,PD0SRMC 1) ) 


EQUIVALENCE 

equivalence 

EQUIVALENCE 


1 

2 

3 

4 

5 

6 
7 


(FOI 601 ) ,AVGNRG ID) 
I EO 1 1 201 ) .NUFNBD I 1) ) 
I E0( 1801 ) ,SIGNFL(1)) 
I EO I 2701 ) .SOTEENID) 
I GO I 295 1 ) ,FVNDRMtl) ) 
I FO I 3301 ) ,UC SFCI 1 )) , 
I FO I 4201 ) .PROTS till® 


,CNS) , I ANS I 3 ) ,FNS ) 
FNRGYN»CASANSI2) ) 

I EO I 301 ) ,EI BNDS 1 1 ) ) , 
I EO I 901 ) ,NFUTXS 1 1 ) ) , 
I EO 1 1 SOI ) ,DT ST 1 1 ) ) , 
(E0I2101) ,UCNFUT(1) ) , 
I EO I 2901 ) ,ND0SRM(1) 1 , 
I FO I 3001 ) ,urPRlMI 1) ) , 
I FOI 7001 ) ,NUFNRG(1 ) ) , 
frnr«;7ni i .onncp- m i i 


1 WRITE (6,2) 

2 FORMAT I lHl ) 

READ (5*3) NOCDS 

3 FORMAT! 12) 

DO 5 M-l, NOCDS 

READ (5.4) (CARD! J) ,J*1 ,12) 

4 FORMAT 1 1 2A6 ) 
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5- WRITE (6,4) (CARD! J) ,J=1,12) 

MOVE = 1 

kntrp=i calculate primary protons only 

KNTRP= 2 CALCULATE PRIMARY PROTONS PLUS FIRST GENERATION SECONDARY 

PROTONS AND EVAPORATION NEUTRONS. 

KNTRP=3 CALCULATE PRIMARY PROTONS PLUS ALL GENERATIONS nF SECONDARY 

PROTONS AND EVAPORATION NFUTRONS. 

KNTRN= 1 CALCULATE FIRST GENERATION CASCADE NFUTRONS 

KNTRN=2 CALCULATE ALL GENERATION CASCADE NEUTRONS AND EVAPORATION 

NEUTRONS. 

read (5,6) spbnd»snbnd,pdsbnd,ndsbnd,bndlow,kntrp,kntrn,sOfeno 

6 FORMAT ( 5F6. 0,314) 

READ (5,6664) KOSW 

6664 FORMAT (3611) 

PROFTD=l .600F-8 

IF (KOSW13 .FQ. 1) PROFTD=PR0FTD*3600. 

INITIALIZE SUBROUTINES AND READ-IN DATA TABLES. 

CALL PROPTYf 1 ,S0FEN0) 

EVNFDO-SUBROUTINE- CALCULATES EVAPORATION NEUTRON DOSF 

CALL EVNFDO (1, DUMMY, DUMMY, DUMMY, DUMMY) 

WRITE (6,6665)BNDL0W,SPBND,SNBND,PDSBND,NDSBND 

6665 FORMAT ( 9X , 55HENERGV CUT-OFF LEVEL FOR INITIAL INCIDENT PROTON FLUX 
A =F7.2,5H MFV./9X.80HFNFRC-Y CUT-OFF LPVFL FOR CALCULATING SFCONDAR 

by particles from incidfnt protons =ft.2.5h mfv./9x,81henergy cut-o 

CFF LEVEL FDR CALCULATING SECONDARY PARTICLES FROM INCIDFNT NEUTRON 
DS =F7. 2 » 5H MfV./9X,56HFNfRgY CUT-OFF LFVFL FOR PROTON DOSF DUE TO 
FI ON I Z AT I ON = F7 . 2 , 5h MEV./9X ,52HENERGY CUT-OFF LFVFL FOR NEUTRON DO 
FSF CALCULATIONS =F7.2,5H MEV. ) 

DOSEK-SUBROUTINE- CALCULATES MASS STOPPING POWFR.RBE, AND FLUX 

TO DOSE CONVERSION FACTOR 

YIELDS-SUBROUTINE-CALCULATES YIELDS OF SFCONDARY PARTICLES FOR 

CASCADF PROTDNS, CASCADE NEUTRONS, AND EVAPORATION NFUTRDNS 

RANGE-SUBROUTINF-CAlCulAtES range from energy or energy from range 

XS-SUBROUTINF-CALCULATFS PROTON AND NFUTRON CROSS-SECT 1 0N$ 

A SNRG- SUBROUT INE-CALCULATES ENERGY OF CASCADE PARTICLES 

INVALU-SUBROUTINE-CALCULATES INCIDENT ENERGY FROM EXIT ENERGY 

• %***I NPlJT AND ALSO CALCULATES THE INPUT PROTON SPECTRUM 
CALL INVALU 
GO TO (11,7) »K0SW7 

.....NONUBD ~ NUMBFR OF NEUTRON FNFRGY BOUNDER IES 
NUENBD- NEUTRON ENERGY BOUND 

7 READ (5*8) NONUBD , ( NUENBD (K ) *K=1 .NONUBD ) 

8 FORMAT ( I3/(7F10.0) ) 

9 DO 10 K = 2, NONUBD 

SIGCNP( K— 1 ,1 )— CASCADF NEUTRONS IN GROUP K-l, FIRST GENERATION 

SIGCNP ( K-l ,1 )=0.0 

SIGCNP(K-l,2)-CASCADE NEUTRONS IN GROUP K-l , SECOND AND HIGHER 

GENERATIONS 

SIGCNP ( K-l, 2)=0.0 

C.....NUENRG- NEUTRON ENERGY(MID-POINT OF INTERVAL) 

10 NUFNRG(K-1)*0.5*( NUENBD ( K )+NUENBD ( K-l ) ) 

GO TO 13 

11 N0NUBD*MAX+1 

DO 12 K=l, NONUBD 

12 NUENBD(K)*FI (K) 

GO TO 9 

13 IPFLTO-0. 

IPDRAD-O. 

!PDREM=0. 

DO 14 J*1,MAX 

C PHI JI { J,1 ) -PRIMARY PROTONS IN ENERGY GROUP J 
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PHI J I ( J » 1 ) = OPPPM(J) 

C PHI JI ( J » 2 ) -FIRST GFNFRATI ON CASCADE PROTONS IN FNERGY GROUP J 

PHI JI ( J»2)*0,0 

PHIJKJ.3) -SECOND AND HIGHER GENERATION CASCADE PROTONS IN GROuPJ 

PHIJI ( J » 3 ) =0 » 0 

EIBND-EX I T FNERGY BOUNDS 

EIBNDS ( J ) «=FT ( J) 

C AVGNRG ( J ) — AVFRAGF FNFRGY OF A GROUP 

AVGNRG ( J)«FTBARt J) 

I PFLTO* I PFLTO+PH I JI ( J,1 ) 

CALL DOSFK ( AVGNRG < J ) ,DEDX »RAB,TFF ,2 ) 

DUMMY = PHIJI (J»l J *DFDX 
IPDRAD * IPDRAD+DUMMY 

14 IPDREM = IPDREM+DUMMY*RABIEF 
IPDRAD * IPDRAO*PR0FTD 
IPDREM = IPDREM*PROFTD 

GO TO (1414,83) ,KOSW5 
1414 MAXP1=MAX+1 

EIBNDS(MAXP1 )*FI (MAXP1 ) 

XX=O,0 

15 N0LAY=0 
SUMENT=0, 

nox- number of large delta x 

DO 82 M*1 ,NOX 
FVNUDS(M) =0,0 

fvndrm(M) » n,o 

IF (M . EO, 1) GO TO 16 

IF (PROPNO(M) ,E0. PROPNO(M-l) )G0 TO 1819 

16 CALL PROPTY ( 3 *M ) 

17 DO 1717 J=2»N0NUBD 

1717 CALL XS ( NUFNRG ( J-l ) ,NEUTXS ( J— 1 ) ,2) 

DO 18 J=1,MAXP1 

IF (FIBNDS(J) ,FQ. 0.0) GO TO 1818 

18 CALL RANGF (FIBNDS(J) ,DIST (J) ,2) 

GO TO 1819 

1818 D I ST ( J ) =0 .0 

C LIMIT- NUMBFR OF SMALL DELTA X 

1819 LIMIT*N0D2X(M) 

HAFD2X- DELTA X/2 

c D 2 x (M )— small delta x 

HAFD2X«D2X(M)*0.5 

19 DO 50 N-l, LIMIT 

C XX-DEPTH INTO THE SHIFLD TO WHICH THE CALCULATION HAS PROGRESSED 

XX*XX+D2X ( M 1 

NOLAY- NUMBFR OF SMALL DELTA X 

NOLAY-NOLAY +1 
XMID( NOLAY)* XX-HAFD2X 
CNTOTS ( NOLAY , 1 ) * 0,0 
CNTOTS ( NOLAY , 2 ) * 0,0 
SOTEENf NOLAY)*0,0 

C ENTOTS(NOLAY)-SUM OF EVAPORATION NEUTRONS PRODUCED IN D2X(M), 

ENTOTS ( NOLAY ) * 0,0 
DO 20 J* 1 »MAX 

UCPRIM ( J ) -UNCOLL IDED PRIMARY PROTONS 

UCPRIM ( J ) *0, 0 
DO 20 JJ"1 , 2 

C UCSFC( J * JJ ) - uncollidfd sfcondary protons 

UCSFC ( J » JJ ) *0, 0 

CPP- CASCADE PROTONS PRODUCED 

20 CPP(J,JJ)«0,0 
DO 23 JJ*1 ,2 
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DO 23 J*2»NONURD 
UCNFUT ( J— 1 * JJ ) *0. 

C CNP ( J— 1 » JJ)- CASCADE NEUTRONS PRODUCFn 

23 CNP ( J-l » JJ ) *0. 0 

C DT ST- PROTON RANGF 

24 DIST(1)=DIST(1)-D2X(M) 

CALL RANGF(DIST(1 ) .FIBNDSd ) *3) 

CALL XStAVGNRGtl ) »SIGNFL(1 ) ,1) 

DO 27 J=2,MAXP1 
DIST( J)»DIST( J)-D2X(M) 

IF (DIST(J) .GT. O.O) GO TO 23 
FIBNDS f J) *0.0 

IF ( ( DT ST ( J-l ) +DI ST ( J ) ) ,GE. 0.0) GO TO 26 
MIDNRGf J— 1)*0.0 
GO TO 2626 

25 CALL RANGE (DIST(J) *FTBNDS(J)*3) 

26 MIDNRGt J-l ) =0.5* (FIBNDS ( J)+E1BNDS (J-l ) ) 

2626 IF(AVGNRO(J) .LF. l.F-5) GO TO 28 

27 CALL XS (AVGNRG(J) ,SIGNFL ( J ) , 1 1 
MIN*MAX 

GO TO 29 

28 MIN * J-l 

29 GO TO (32*30) ,K0SW9 

30 WRITE (6,31) (AVGNRG(NN) .SIGNEL(NN) ,NN»1,MIN) 

31 FORMAT (1HL8X.15HCR0SS— SECTIONS, /1HJ12X*6HENFRGY13X»9HINELASTIC/ 

1 (2X,1P2E20.5) ) 

32 DO 3838 J=1,MIN 

C DUMMY-THIS NAME IS USED FOR SEVERAL DIFFERENT QUANTITIES IN THE 

C CALCULATION of SECONDARY particle production, particle attenuation, 

c and particle dose calculations 

CALL RANGE ( AVGNRG ( J ) , DUMMY , 2 ) 

IF (DUMMY .LE. D2X(M)) GO TO 34 
DEGMLT=EXP (-S IGNEL ( J) *D2X ( M ) ) 

UCPRTM ( J ) *PHI JT ( J»1)*DFGMLT 
IF (KNTRP .EG. 1) GO TO 35 
DO 33 JJ»2,KNTRP 

33 UCSEC ( J, JJ— 1 ) »PHI JI ( J » JJ) *DEGMLT 
GO TO 35 

34 UCPRIM{ J) *0.0 

IF (DUMMY .LT. HAFD2X) GO TO 3838 

C DEGMLT- PROTON ATTENUATION FACTOR ACROSS DELTA X OR DELTA X/2 

DEGMLT«EXP(-SIGNEL( J)*HAFD2X) 

35 IFdKNTRP .FQ. 1 .AND. KNTRN .FQ. 0) .OR. (AVGNRG (J) ,LT. SPBND) ) 
AGO TO 3838 

C COLFRA- CDLLIDFD FRACTION 

36 COLFRA » 1,0 - DEGMLT 
DO 3738 N0M*1 ,NOFCOM 

C ANS(K) - YIELDS OF SECONDARY PARTICLES - K-l , CASCADE PROTON — 

C K-2, CASCADE NFUTRON — K«3, EVAPORATION NEUTRON. 

CALL YIELDS ( AVGNRG ( J ) ,ANS ,NOM , 2 ) 

CALL CASNRG ( AVGNRG ( J ) ,CASANS ,NOM , 2 ) 

IF ( KNTRP .LT. 2) GO TO 3737 

C ENRGYP- ENERGY OF CASCADE PROTON 

C DUMRNG- RANGE OF SECONDARY PROTON AT BIRTH 

CALL RANGE (ENRGYP, DUMRNG, 2) 

DUMRNG-DUMRNG-HAFD2X 
IF (DUMRNG .LE. 0.0) GO TO 37 
CALL RANGF ( DUMRNG , PROF ,3 ) 

CALL XS (FNRGYP,ATXSFC,1) 

DUMMY«C0LFRA»CPS*EXP(-ATXSEC*HAFD2X) 

CALL SORT (PROE ,EIBNDS,MAX ,IND1 
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CPP( IND*1 )*PHI JI < J»1 )*DUMMY +CPP ( IND.l ) 

SUMaPHI JI ( J,1 ) 

IF (KNTRP .LF. 2) GO TO 37 

CPP( IND» 2) *DUMMY* (PHI JI ( J.2J+PHTJI (J,3) > +CPP ( IND,2) 

SUM=SUM+PHIJT ( J » 2 )+PHT JI (J»3) 

37 IF ( (PROPNO(M) .LT.50) .AND. (NOM.FQ.l )) GO TO 3738 

ENTOTS(NOLAY)-EvAPORATION neutrons produced by all protons 

ENTOTS(NDLAY) = ENTOTS(NOLAY)+COLFRA*ENS*Sl)M 

3737 IF ( KNTRN .EQ. 0) GO TO 3738 

FNRGYN- FNERGY OF CASCADF NFUTRON 

CALL XS (FNRGYN »ATXSEC ,2 ) 

ATTFN = FXP ( -ATXSEC*HAFD2X ) 

DUMMY = PHTJI ( J,1 )*COLFRA*CNS 

CALL SORT ( FNRGYN, NUENBD.NONUBD-1 ,IND) 

CNT0TS(N0LAY,1 ) = CNT0TS(N0LAY»1 I+DUMMY 

CNP ( I ND * 1 ) -CASCADE NEUTRONS PRODUCED BY PRIMARY PROTONS 

CNP ( I ND » 1 ) * CNP( IND.l ) +DUMMY*ATTEN 
IF (KNTRN .EQ. 1) GO TO 3738 

C CNP( IND.2) -CASCADE NEUTRONS PRODUCED BY SECONDARY PROTONS 

DUMMY = (PHI JI ( J,2)+PHTJl(J»3n*C0LFRA*CNS 
CNT0TS(N0LAY,2) = CNTOTS (NOLAY ,2 I+DUMMY 
CNP ( I ND * 2 ) * CNP ( I ND * 2 I +DUMMY*ATTFN 

3738 CONTINUE 
3838 CONTINUE 

IF (KNTRN .FQ. 0 .OR. XX .EQ. D2X(1)) GO TO 43 
DO 42 J*2 .NONUBD 

DUMMY= C XP (-NFUTXS ( J-l ) *D2X(M)) 

DO 39 JJ=1, KNTRN 

39 UCNFUT ( J— 1 » J J ) =SIGCNP( J-l , JJ)*DUMMY 

IF ((KNTRN .LT. 2 .AND. KNTRP .LT. 3) .OR. (NUENRG(J-l) .LT. SNBND 
A) ) GO TO 42 

SUM=SIGCNP(J-1,1 )+S IGCNP ( J-l *2 ) 

I F ( SUM .FQ. 0.0) GO TO 42 

C0LFRA=1. -DUMMY 

DO 4141 N0M=1»NPFC0M 

CALL YIELDS (NUFNRG( J-l ) ,ANS,N0M,3 ) 

CALL CASNRG (NUFNRG( J-l ) .CASANS »NOM , 3 ) 

IF (KNTRP .LT. 3) GO TO 40 
CALL RANGE ( FNRGYP»0UMRNG»2 ) 

DUMRNG=DUMRNG-HAFD2X 
IF (DUMRNG .LE. 0.0) GO TO 40 
CALL RANGE (DUMRNG, PR0E.3) 

CALL XS ( FNRGYP ,ATXSEC»1 ) 

DUMMY=COLFPA*CPS*EXP ( -ATXSFC*HAFD2X ) 

CALL SORT (PROF ,FIBNDS,MAX ,IND) 

CPP ( IND*2 )»DUMMY*SUM +CPD(IND,2) 

40 IF (KNTRN .LT. 2) GO TO 4141 

CALL XS (FNRGYN, ATXSFC, 2) 

DUMMY = SUM*COLFRA*CNS 

CALL SORT (FNRGYN, NUENBD, NONUBD-1 ,IND) 

C CNP( IND»2)-CASGADF NFUTRONS PRODUCED RY CASCADF NFUTRONS 

CNT0TS(N0LAY,2) * CNT0TS(N0LAY,2 I+DUMMY 

CNP ( I ND ,2 ) = CNP ( IND.2 )+DUMMY»FXP(-ATXSEC*HAFD2X) 

ENTOTS(NOLAY)-FVAPORATlON NEUTRONS PRODUCED BY CASCADE NEUTRONS 

41 ENTOTS(NOLAY) a ENTOTS(NOLAY)+COLFRA*FNS»SUM 
4141 CONTINUE 

42 CONTINUE 

43 DO 45 J=1,MAX 

PHI JI ( J,1 )*UCPRIM( J) 

IF (KNTRP .FQ. 1) GO TO 45 
DO 44 JJ*2, KNTRP 


E-3080 


« 


51 


44. D HI JI ( J.JJ)*UCSFC( )+CPP ( J.JJ-l ) 

45 AVGNRGt J)*=MIDNRG( J) 

IF (XNTRN .FQ. 0) GO TO 47 
DO 46 JJ*1 , KNTRN 
DO 46 J=2,NnNUBD 

46 S IGCNP ( J-l . JJ ) =L)CNFUT ( J— 1 ,JJ)+CNP( J-l ,JJ) 

GO TO 48 

47 IF (XNTRP .FQ. 1) GO TO 50 

48 SUMFNT=SUMFNT+FNTOTS( NOLAY) 

CNTOTS ( NOLAY * 1 ) «CNTOTS (NOLAY » 1 ) /D2X ( M ) 
CNTOTS(NOLAY,2)=CNTOTS(NOLAY,2) /D2X(M) 

SOTFFN ( NOLAY ) =FNTOTS ( NOL AY ) /D 2X ( M ) 

CALL FVNFDO ( 2, NOLAY. M,N, DUMMY) 

50 CONTINUE 

DO 51 X = 1 .9 

51 TOTALS(M,X)«0.0 

SUMMATION OF PROTON PARTICLES. 

52 DO 53 X=l, XNTRP 
DO 53 MM=1.M1N 

TOTALS (M.X) « T0TALS(M,X)+PHUI (MM.X) 

IF (X.LT.2) GO TO 33 

TOT ALS ( M , X+2 ) * T0TALS(M,X+2 )+CPP (MM.X-1 ) 

53 CONTINUE 

54 IF (XNTRN .FQ. 0) GO TO 56 

C SUMMATION OF NFUTRON PARTICLFS. 

DO 55 K=l. KNTRN 
DO 55 MM-2.N0NUBD 

TOTALS (M.X+5) * TOTALS (M.X+5 ) +S IGCNP (MM-1 ,X ) 

55 TOTALS (M.X+71 * TOTALS (M.X+7 ) +CNP ( MM-1 ,X ) 

C.... .PROTON DOSF CALCULAT I ONS. 

56 DO 57 XK-1,3 
DO 57 X*1 .5 

C PDOSF ( M *X »XX ) - PROTON DOSF IN RAD 

PDOSE (M.X. XX ) * 0,0 
C.....PDOSRM(M,X,XX)-PROTON DOSE IN RFM 

57 PDOSRM(M.X.XX) a 0.0 
DO 59 MM*1 .MIN 

C DFDX-MASS STOPPING POwfR OF PROTONS 

C RABIEF-RELATIVE BIOLOGICAL EFFECTIVENESS ( RBE ) 

CALL DOSEK ( AVGNRG(MM) .DEDX.RABIEF .2) 

DO 59 K«1»KNTRP 
DUMMY » PROFTD*PROTS(MM,X)*DEDX 
PDOSE (M.X.l) - PD0SE(M.X.1)+DUMMY 
PD0SRM(M»X*1 ) > PDOSRMtM.X.l ) + DUMMY*RABIFF 
GO TO (58. 5757) .X0SW21 
5757 IF (AVGNRG(MM) .LT. 10.) GO TO 5ft 

XOFF- NUCLFAR DOSF VARIABLE 

CALL DOSFK (AVGNRG(MM) ,X0FE*DUMMY.4) 

DUMMY - XOFF*PROTS(MM.X> 

PDOSF (M.X. 2) » PDOSF ( M »X *2 ) + DUMMY 
PDOSRM ( M. X, 2) * PDOSRM (M » X ,2 ) + DUMMY*RABIEF 
50 IF (X.LT.2) GO TO 59 

DUMMY * PROFTD*PROTS(MM,X+2)*DFDX 
PDOSF (M.X+2.1) * PDOSF ( M ,K+2 . 1 ) + DUMMY 
PD0SRM(M,X+2.1) = PDOSRM (M.X+2.1 ) + DUMMY*RAPTFF 
GO TO (59.5858) .X0SW21 
5ft5ft I F ( AVGNRG ( MM ) .LT. 10.) GO TO 59 
DUMMY * X0FF*PR0TS(MM,X+2) 

PDOSF (M.X+2.2) * PDOSF (M.X+2 .2 ) + DUMMY 
PDOSRM ( M ,X+2 . 2 ) = PDOSRM (M.X+2 . 2 ) + DUMMY*RAB IFF 
59 CONTINUE 
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GO TO < 5960,5958 ) ,KOSW21 

5958 DO 5959 <*1*5 

PDOSF (M,K,3) = PDOSF (M,K,1) + PDOSF(M,K,2) 

5959 PD0SRM(M,K,3) = PDOSRM(M,K,l ) + PDOSRM (M ,K . 2 ) 

C,. # . .NEUTRON DOSF CALCULATIONS. 

5960 DO 60 KM .4 

C.....NDOSF(M*KK)- DOSF AT THICKNESS X FROM CASCADE NFUTRONS IN RAD 
NDOSF (M.K ) = 0.0 

C • • . . .NDOSRM ( M.KK ) -DOSE AT THICKNESS X FROM CASCADF NFUTRONS IN REM 

60 NDOSRM f M * K ) = 0.0 

IF (KNTRN.FG.O) GO TO 6161 
DO 61 K»2 *NONUBD 

CALL DOSFKfNUENRG(K-l) ,CNK1,CNK2,3) 

DO 61 KKM.KNTRN 

NDOSF (M.KK) = NDOSF ( M , KK )+NFWTS ( K-l ,KK ) *CNK 1 
NDOSE ( M.KK+2 ) * ND0SF(M,KK+2)+NEWTS(K-l.KK+2)*CNKl 
NDOSRM (M.KK) * ND0SRM(M,KK)+NFWTS(K-1 ,KK)*CNK2 

61 ND0SRM(M,KK+2) = NDOSRM ( M ,KK+2 ) +NEWTS ( K— 1 »KK+2 ) *CNK2 
GO TO 62 

6161 IF (KNTRP .FQ. 1) GO TO 6262 

62 CALL FVNFDO ( 3 ,NOLAY ,M,L I MIT, FVNUDS ( M ) ) 

FVNDRM(M) = 10.0*FVNUDS(M) 

6262 GO TO ( 82 ,63 ) ,K0SW11 

C OUTPUT OF ALL DATA. 

C 

C.... .PRIMARY AND SFCONDARY PROTON 

63 KNTPG = 0 
ASSIGN 64 TO LIMF 

64 LINF1* KNTPG*54+1 
KNTPG = KNTPG +1 
LASTLN = KNTPG*54 

IF (LASTLN - MAX ) 67,66,65 

65 LASTLN * MAX 

66 ASSIGN 69 TO LIME 

67 WRITE (6,681 XX 

68 FORMAT! lH 146 X, 21 HSH1 ELD THICKNFSS , X *F8.2,9H GM/CM**2/1H060X , 34HP 
1R0T0N FLUX SPECTRA— PR0T0NS/CM**2 ) 

IF ( KOSW ( 1 3 ) .FQ. 1 )WRITF(6,6868 ) 

6868 FORMAT ( 1H+94X ,4H-SFC ) 

WRITf (6,6869) ( J , AVGNRG ( J) ,(PROTS(J,K) ,K = 1,5) , J='. I NFl , LASTLN ) 

6869 FORMAT ( 12X,6HFNFRGY13X,7HPRIMARY13X,34H - SFCONDARY PROTONS AT 

IX 9X,35H SFCONDARY PROTONS IN DX /13X ,3HMFV15X,7HPR 

20T0NS3X.2 ( 1 IX , 1 OHF I R ST GEN. 12X , 1 1 HHIGHFR GFN. ) / ( 17 ,OPF12.3 , 1P5E22. 
36) ) 

GO TO LIMF, (64,69,74,79) 

69 WRITF (6,70) ( TOT ALS ( M ,K ) ,K=1 , 5 ) , ( PDOSF (M ,K , 1 ) ,K-1 , 5 ) 

70 FORMAT ( 1HJ3X,6HT0TALS9X,1P5F22.6/1HJ3X,9HD0SF — RAD6X .5F22.6 ) 

IF ( KOSW (13) .FQ. 1) WRITF (6,71) 

71 FORMAT ( lH+1 2X » 4H/HP ) 

WRITE (6,72) . (PD0SRM(M,K,1),K«1,5) 

72 FORMAT (4X,9HD0SF — RFM6X , 1 P5F22 .6 ) 

IF ( KOSW( 13 ) .FQ. 1) WRITE (6,71) 

IF (KNTRN .FQ. 0) GO TO 82 

C CASCADE NEUTRON OUTPUT. 

73 KNTPG « 0 
ASSIGN 74 TO LIMF 

74 LI NFl = KNTPG*54+1 
KNTPG = KNTPG + 1 
LASTLN ■ KNTPG#54 

IF (LASTLN - NONUBD+1 ) 77,76,75 

75 LASTLN - NONUBD-1 
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76- ASSIGN 79 TO LIME 

77 WRITE (6,78) XX 

78 FORMAT ( 1H133X.21H SHIELD THICKNESS * X =F8.2,9H FM/CM**2/1H042X,44HC 

iascaof neutron flux spectra — nfutrons/cm**2) 

IF ( KOSW ( 13 ) .FO. 1 )WRITF(6,7878 > 

7878 FORMAT ( 1H+86X »4H _ SFC ) 

WRITF (6,7879) ( J.NUFNRgU) , (NFWTS(J,K) ,K = 1,4),J = LINF1,LASTLN) 

7879 F0RMAT(12X,6HENFRGY11X,33H- CASCADF NFUTRONS AT X - 11X.34H 

A - CASCADE NFUTRONS IN DX -/I 3X ,3HMFV3X , 2 (11X , 10HF I RST GEN. 

B12X.11 HH I GHE R GEN . ) / ( 17 , OPF1 2 .3 , 1 P4E22 .6 ) ) 

GO TO LIME, (64,69,74,79) 

79 WRITF (6,80) ( TOTALS ( M , K ) ,K=6 ,9 ) , ( NDOSF ( M ,K ) , K«1 ,4 ) 

80 FORMAT ( 1HJ3X,6HT0TALS9X,1P4E22.6/1HJ3X,9HD0SF— RAD6X.4E22.6) 

IF (K0SW(13)- .FO. 1) WRITE (6,71) 

WRITF (6,72) (NDOSRM(M.K) ,K=1,4) 

IF ( KOSW ( 13 ) .FO. 1) WRITF (6,71) 

82 CONTINUE 

C DOSE TABLES — RAD OR RAD/HR 

83 M = 0 

LINKNT=26/K0SW21 
KONGIB * 2#KOSW21-l 

84 WRITF (6,85) 

85 FORMAT ( 1H159X .9HD05F — RAD) 

I F ( KOSW( 1 3 ) . FQ. 1) WRITE (6,8585) 

8585 FORMAT ( 1H+68X »3H/HR ) 

WRITE (6,86) (LLL,LLL=1 ,6) , I PDRAD , IPDRAD , IPDR AD 

86 FORMAT (77H0SHIFLD PRIMARY — SECONDARY PROTON CASCADE N 

1FUTR0N EVAPORATION ,5 ( 5HT0TAL6X ) /17H THICKNESS PROTON , 2 ( 5X , 5HF IR 

2ST6X .6HHIGHFR) ,5X .59HNEUTR0N SECONDARY CASCADE PROTON N 

3FUTR0N D0SE/9H GM/CM»*21 3x ,4 ( 1 1HGFNFRATI0N ) ,11x,6HPR0t0N5x,7HN 
4EUTR0N/12X,6(3H ( 1 1 , 1H ) 6X ) , 53H( 2 ) + ( 3 ) (4)+(5) ( 1 ) + ( 2 > + ( 3 ) ( 4 ) + 

5 ( 5 ) + ( 6 ) (1)THRU(6)/8H0 0.00,1PE13.4,F88.4,F22.4) 

GO TO (87,90) ,K05W5 

87 M = M+l 

TCN = NDOSF t M , 1 ) + ND0SE(M,2) 

TN * TCN+FVNUDS(M) 

DO 88 KK-1, KONGIB 

TSP(KK) « PD0SF(M,2,KK) +PDOSF (M» 3 »KK ) 

TP(KK) * TSP(KK)+PD0SF(M,1,KK) 

88 TD(KK) = TP ( KK ) +TN 

WRITE (6,89) X (M ) , ( PDOSE (M ,K , 1 ) ,K«1»3) , ( NDOSF ( M ,K ) ,K*1 ,2 ) .FVNUDS 
1 ( M) , TSP ( 1 ) ,TCN,TP(1) ,TN,TD(1 ) 

89 FORMAT ( 1H0F7.2*2X*1P11f11.4) 

GO TO (8991,8989) ,K0SW21 

8989 WRITE (6,8990) ( ( PDOSE ( M ,K , N ) ,K-1 , 3 ) ♦ TSP ( N ) , TP ( N ) ,TD ( N ) ,N*2 ,3 ) 

8990 F0RMAT(9X,1P3E11.4,F44.4,2F22.4) 

8991 IF (M.GE.NOX) GO TO 90 

IF (MOD(M.LINKNT)) 87,84,87 
C DOSE TABLES — REM OR REM/HR 

90 M * 0 

91 WRITF (6,9?) 

92 FORMAT ( 1H159X ,9HD0SF— REM) 

IF (K0SW13 .FO. 1) WRITF ( 6,8585 ) 

WRITE (6,86) (LLL,LLL«1,6),IPDREM,IPDREM,IPDREM 
GO TO (93,95) .K0SW5 

93 M « M+l 

TCN ■ ND0SRM(M,1 )+ND0sRM(M,2 ) 

TN » TCN+EVNDRM(M) 

DO 94 KK*1, KONGIB 

TSP(KK) » PDOSRM(M,2»KK)+PDOSRM(M,3,KK) 

TPtKK) = TSP ( KK ) +PDOSRM ( M, 1 , KK) 
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94 TD(KK) * TP { KK ) +TN 

WRITE (6,89) X(M) ,(PD0SRM(M,K,1 ) ,K»1,3) ,(NDOsRM(M,K) ,K» 1 ,2) »EVNDRM 
1 (M) ,TSP(1) ,TCN,TP(1 ) *TN ,TD ( 1 ) 

GO TO ( 9495 *9494 ) »KOSW21 

9494 WRITE (6,8990) ( ( PDOSRM (M ,K ,N ) ,K*1 ,3 ) ,TSP ( N ) , TP ( N ) ,TD ( N ) ,N*2 ,3 ) 

9495 IE (M.GF.NOX) GO TO 95 

IF (MOD(M,LINKNT) ) 93,91,93 

95 GO TO (100,96) »K0SW17 

: FLUX TABLFS 

96 WRITE (6,9696) 

9696 FORMAT ( 1H148X »21HFLUX— — PARTICLFS/CM**2 ) 

IF (K0SW13 ,F0, DWRITF (6,9697) 

9697 FORMAT ( 1H+69X ,4H-SFC ) 

WRITE (6,9698) (LLL ,LLL = 1 ,5 > , IPFLTO 

9698 FORMAT ( 7H0SHIELD8X »7HPRI MARY8X ,6 OH SECONDARY PROTON 

1 CASCADE NEUTRON , 2 ( 5X , 5HT0TAL6X ) / 10H THICKNES C 5X ,6HPR0 

2T0N10X,2( 5HFI RST1 IX ,6HHIGHER10X ) ,9HSEC0NDARY7X,7HCASCADE/9H GM/CM* 
3*222X»4(1 OHGENERATI 0N6X ) ,6HPR0T0N10X ,7HNEUTR0N/17X , 5 ( 3H ( I1,1H)1 
41 X) ,7H(2) + (3)9X,7H(4) + (5) /1H03X ,4H0.001PF15.5 ) 

GO TO (9797,1) *K0SW5 
DO 98 M=1,N0X 

TSP * T0TALS(M,2)+T0TALS(M,3) 

TCN = T0TALS(M,6)+T0TALS(M,7) 

WRITE (6,99) X ( M ) ♦ (TOTALS (M,K) ,K»1 ,3) , ( TOTALS < M ,K ) ,K=6»7) ,TSP(1 ) , 


9797 


98 


1TCN 

99 FORMAT ( 1H0F7,2,2X»1P7F16.5) 
ion go in d,ioi) ,r.nswi 
C NEUTRnN SOURCE TERMS TABLF 

101 LETA =N0LAY/2 
IDXSC = (NOLAY+1 )/2 
WRITF (6,102) 

102 FORMAT ( 1H1 ,2 ( 9X ,6HSHIELD10X ,32HNEUTRON SOURCE TERM — NEUTR0NS/GM7X ) 
A) 

IF(K0SWl3 ,F0. 1) WRITF(6, 10202 ) 

10202 FORMAT ( 1H+2 ( 57X .AH-S^C , 3X ) ) 

WRITE (6,10203) 

10203 FORMAT ( 2X * 2 ( 6X »9HTHICKNES55X* 11HFVAP0RAT I0N6X ,27H - CASCADE 

A ) /3X,2(6X,10H(GM/CM**2)21X,10HFIRST GEN ,6X , 1 1HHIGHFR GFN 


B. ) ) 


DO 103 K»1,LETA 
NDX= IDXSC+K 


103 WRITE (6,104) XMID(K) ,SOTEFN(K) , (CNTOTS(K.J) ,J*1,2) ,XMID(NDX) , 
ASOTFFN(NDX) , (CNTOTS(NDX,J) , J»1 , 2 ) 

104 F0RMAT(F16.3,1P3F17.5,0PF13.3,1P3E17.5) 

I F ( MOD ( NOLAY , 2 ) ,EQ. 1) WRI TE ( 6 , 1 04 ) XMID (LFTA+1 ) ,SOTEFN ( LETA+1 ) , 
A(CNT0TS(LETA+1.J) ,J=1,2) 

GO TO 1 
END 

SIBFTC EVNFDO LIST ,DFCIC 


SUBROUTINE FVNFDO ( INDFX 
EVAPORATION NEUTRON DOSE 

,layno,mm,nn,dose> 

calculation 

COMMON 

D2X ( 20 ) 

♦ 

MAX , 

EO ( 300 , 20 ) 

COMMON 

FT (300) 

* 

DFK300) , 

FIRAR (300) 

COMMON 

OP ( 300 ) 

9 

0PPRM(3OO) , 

N0D2X ( 20 ) 

COMMON 

X( 20) 

9 

NOX , 

FNTOT S ( 200 ) 

COMMON 

DX ( 20 ) 

9 

PROPNO ( 20 ) , 

Cl 

COMMON 

COMMON 

PDSBND 

9 

ndsbnd , 

ENRG(IOO) , 

BNDLOW 

RNG(IOO) 

COMMON 

EBOMBPI 25*4 

) 9 

FBOMBNt 25 ,4) , 

CPSP( 25,4) 

COMMON 

CPSN (25,4) 

9 

CNSP (25,4) , 

CNSN (25,4) 

COMMON 

ENSP (25,4) 

9 

ENSN (25,4) , 

EBOMP (25,4) 


0010 
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COMMON 

EBOMN (25,4) 

, EPROP (25,4) 

, EPRON (25,4) 

COMMON 

ENEUP (25,4) 

, ENEUN (25,4) 

, ENRGP ( 25 ) 

COMMON 

ENRGN ( 75 ) 

, XSMBP ( 25 ) 

, XSMBN ( 75 J 

COMMON 

SNRG(IOO) 

, RBENRG ( 20 ) 

, C1NRG ( 40 ) 

COMMON 

C2NRG ( 40 ) 

, CNRG( 2 ) 

, SOFF (100) 

COMMON 

RBF( 20) 

, C0NK1 (40) 

, CONK? ( 40 ) 

COMMON 

CONK ( 2 ) 

• LFNGTH(8) 

, GMWT 

COMMON 

LSOFE 

, LRBE 

, LK1 

COMMON 

LK2 

, LK 

, NOFCOM 

COMMON 

MOVF 

, KOSW( 36 ) 

* KONSWT 


EQUIVALENCE ( KOSW ( 13 ) .K0SW13 ) 

DIMENSION P ( 20 ) »H ( 20) »BNDANG( 11 ) , RADANG ( 10 ) ,COSCON ( 10 ) .COSRAD ( 10 ) , 
AS (20) ,SUMGPA(10) .SUMSR ( 10,200) ,SUMHRP ( 1 0 , 200 ) ,K(20) 

INTEGER PROPNO,PTEST,PDIE 

REAL KS.KSR 

REAL LI ,KI ,KLSR,K 

DATA C2»C3,C4/0. 3492, 0.4223, 0.6984/ 

GO TO ( 1,20,27) , INDEX 

1 READ (5,3) NOX.NOANG, (X(J) ,N0D2X( J) ,PROPNO( J) ,P( J) ,H(J) ,S( J) ,K( J) , 
AJ=1 ,NOX) 

3 FORMAT(2I4/(F6.0,2I4,3E8.0,F6.0) ) 

WRITE (6,4) 

4 FORMAT ( 1H09X.81HX.MIN X»MAX NUMBER OF MATERIAL DENSITY 

A HYDROGEN RFMOVAL XSECT/12X,77H(GM/CM**2) INCREMENTS 

B NUMBER (GM/CM**3) RATIO (CM**2/GM)) 

XMIN=0.0 
DO 2 J=1 ,NOX 

WRITE (6,5) XMIN,X( J) ,N0D2X(J) ,PROPNO( J) ,P(J) ,H(J) ,S(J) 

5 F0RMAT(F15.2,E9.2,I10,I12,F14.5,F12.3,F15.4) 

2 XMIN=X(J) 

WRITE (6,6) NOANG 

6 FORMATf 1H08X, 59HNUMBFR OF ANGLFS IN EVAPORATION NEUTRON DOSE CALCU 
1L AT IONS =13) 

ENFTD = 5.389F-9 

IF (K0SW13 .EQ. 1) FNFTD=FNFTD*3600. 

ML A ST = 0 

lsthyd=o 

lsthla=o 

C.....DELANG- DELTA ANGLE 


DELANO = 1.5707963/FLOAT ( NOANG ) 0180 

C BNDANG- BOUNDARY ANGLES OF DELTA ANGLE 

BNDANG ( 1 ) = 0.0 0190 

C RADANG- MID-POINT OF ANGLE INTERVAL 

RADANG(l) ■ 0.0 0200 

BNDANG ( 2 ) = DFLANG 0210 

CANGL1 * Cns( BNDANG (2 ) ) 0220 

C COSCON- DELTA SOLID ANGL C / { 4*PI F*COSR AD ) 

COSCON ( 1 ) * (1.0-CANGL1 J/2.0 0230 

C COSRAD- COSINE OF MID-POINT ANGLE 

COSRAD(l) » 1.0 0240 

IF (NOANG .LT. 2) GO TO 1 

DO 10 J»2, NOANG 0250 

BNDANG ( J+l ) « BNDANG ( J ) +DFLANG 0260 

RADANG ( J ) * ( BNDANG ( J ) +BNDANG ( J+l ) )/2,0 0270 

COSRAD ( J ) = COS (RADANG ( J) ) 0280 

CANGL2 * COS(BNDANG( J+l ) ) 0290 

COSCON ( J ) = (CANGL1-CANGL2) /2 .0 /COSRAD (J ) 

10 CANGL1 » CANGL2 0310 

GO TO 100 0320 

20 kount=layno 

M=MM 
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N=NN 

KM1=K0UNT-1 

IF ( M . FQ. MLAST) GO TO 6 
MLAST-M 

I F ( (M .GT. 1) .AND. (PROPNO(M) .EG. PROPNO(M-l))) GO TO 6 
JUMPO el 

IF (PROPNO(M) .GT. 100 ) (50 TO 150 
C.....H(M)-RATI0 OF HYDROGFN IN MATERIAL M TO THAT IN WATER 

C P (M) — DFNSITY OF MATFRIAL M IN gPaMS/CM**3 

HOPeH(M) /P fM ) 

JUMPO = 2 

C S(M)-REMOVAL CROSS section 

150 KS=K(M)*S(M) 

ASSIGN 120 TO LOKSUN 
IF (M . EO. 1) GO TO 54 
KK=M-1 

42 IF (PROPNO(M) .NE. PROPNO(KK)) GO TO 45 
50 IF (KK .FO. 1 ) GO TO 54 
KK=KK-1 
GO TO 42 

45 PTEST s ( PROPNO (M ) +PROPNO ( KK ) ) /2 

I F ( ( PTEST .GT. 100) .OR. (PTEST .LT. 50)) GO TO 54 
PDIF +PROPNO ( M ) -PROPNO ( KK ) 

IF (PDIF) 45,50,52 
4P L0W=LSTHYD+1 
I GH=M- 1 

layknt=o 

IF ( LSTHYD .NF. 0) LAYKNT *LSTHL A 
MIN=LAY<NT+1 
DO 160 INFLOW, IGH 
NUMB=N0D2X( IN) 

D2XKS=02X<IN)*K(IN)*S( IN) 

DO 160 JJ=1,NUMB 
LXMlsLAYKNT 
LAYXNTeLAYXNT+1 
DO 160 N0A*1,N0ANG 
KSR*D2XKS/C0SRAD(N0A) 

S UMSR (NOA, LAYKNT) -SUMSR ( NOA, LAyKNTJ- 0.5*KSR 
IF(LAYKNT .FO. MIN) GO Tn 153 
DO 151 L0 »min,LKM1 

151 SUM SR ( NO A ,L0) eSUMSR ( NO A *L0) -KSR 

153 IF(LSThY0 .FO. 0) GO TO 152 
DO 154 L0*1,LSTHLA 

154 SUMSR (N0a,L0)=SUMSR( NOA ,LD)-KSR 

152 SUM=HUN) /P( IN)*(0.5+FLOAT(NOD2X( INJ-JJ) )* 
AD2X ( IN) /COSRAD(NOA) 

I NP1=I N+l 
DO 155 KA=INP1,M 

155 SUM=SUM+H(KA)/»(KA)*DX(KA)/COSRADtNOA) 

157 L 1*0. 5+SHM/l 5 . 

IFILI .GT. 1,0) LI-1.0 
KLSR=KSR*LI/K ( IN) 

SUMSR ( NOa .LAYKNT ) -SUMSR ( NOA .LAYKNT ) +0 . 5*KLSR 
IF (LAYKNT .LF. 1) GO TO 160 
DO 159 L0=1,LKM1 

159 SUMSR (NOa,LO)=SUMSR( NOA ,LO)+KLSR 

160 CONTINUE 

lsthla=kount--i 

TF { K ( M ) .EG. 1 . ) GO TO 54 
ASSIGN 130 TO LOKSUN 
GO TO 60 
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52- LI=K(M) 

GO TO 60 
5 4 L 1= l.n 
60 DO 40 J=1,N0ANG 

C RI- SLANT PATH LFNOTH ACROSS DFLTA X 

RI=D2X(M) /COSRADt J) 

KSR=KS*RI 

SUMSR ( J .KOtJNT ) =KSR*0. 5 
IF (KOtJNT .FO. 1) GO TO 174 
IF ( PROPNO(M) .LT. 50) GO TO 129 
LSTP1*LSTHLA+1 

IF ( LSTP1 .GT. KM1 ) GO TO 120 
DO 2525 K0=LSTP1.KM1 
2525 SUMSR (J.KO) =SUMSR ( J »K0 )+KSR 
IF(LSTHLA .FQ. 0) GO TO 174 
129 GO TO LOKSUN, (120.130) 

150 S(JM*H0P* t 0,5+FLOAT (N0D2X (M)-N) )*RT 

67 LI=0.5 +SUM/15. 

IF (LT .GT. 1.0) LI=1.0 
120 KLSR=LI*KSR/K(M) 

DO 25 K0=1, LSTHLA 
25 SUMSR (J.KO) =SUMSR (J,KO)+KLSR 

174 GO TO ( 190,175) .JUMPO 

175 HRP=HOP*R I 

SUMHRP ( J.KOUNT ) =HRP*0 . 5 
IF ( (COUNT .FQ. 1) GO TO 4 
DO 177 K0»1.KM1 

1 77 SUMHRP ( J » KO ) = SUMHRP ( J * KO ) 4-HRP 
Gn TO 40 

190 SUMHRP (J.KOUNT) =0.0 
40 CONTINUE 

I F ( PROPNO ( M ) .LT. 50) LSTHLA=KOuNT 

IF{ (PROPNO(M) .GT. 100) .OR. (N .NF. N0D2X(M))) GO TO 100 

LSTHYD=M 

GO TO lOo 

27 DOSE = 0. 

KOUNT=LAVNO 

M=MM 

lstpi«lsthla+i 

DO 28 J-l.NOANG 
SUMA*0 • 0 
SUMB=0 . 0 

IF(LSTHLA .FO. 0) GO Tn 30 
DO 29 K0=1 .LSTHLA 

IF ( SUMHRP ( J ,K0 ) .LT. 2.) GO TO 202 

TERMA * SUMHRP( J,K0)**C2*ExP(-C3*SIJMHRP( J.KO) »*C4) 

GO TO 200 

202 TFRMA * .772 - ,065*SUMHRP ( J.KO ) 

200 TERMB = FXP ( -SUMSR ( J.KO ) ) 

29 SUMA=SUMA+TFRMA*TFRMB*FNTOTS(KO ) 

SUM A= SUM A *C0 SC (IN ( J ) 

IF (LSTHLA .FQ. KOUNT) GO TO 32 

30 DO 31 KOeLSTPl .KOUNT 

TFRMB » .772*FXP(-SUMSR( J.KO) ) 

31 SUMB=SUMB+TFRMB*ENTOTS(KO) 

SUMB»SUMB*COSCON ( J ) 

32 CONTINUE 

28 DOSF = DOSF +FNFTD* ( SUMA+SUMB ) 

100 RETURN 0650 

FND 0660 

SIBFTC XS LIST. DECK 
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SUBROUTINE XS(F»XSFCT,TNDFX ) 

C CROSS-SECTION CALCULATIONS. 

C INPUT CROSS-SECTION DATA IN MILLI -BARNS. 


COMMON 

D2X ( 20 ) 

* 

MAX 

* 

FO ( 300, 20 ) 

COMMON 

FI ( 300 ) 

9 

DFT (300) 

9 

FIBAR (300) 

COMMON 

OP ( 3 DO ) 

9 

OPPRM (300 ) 

♦ 

N0D2X (20) 

COMMON 

X ( 20 ) 

9 

NOX 

9 

ENTOTSI 200 ) 

COMMON 

DX ( 20 ) 

9 

PROPNOf 20) 

9 

Cl 

COMMON 

PDSBND 

9 

ndsbnd 

9 

BNDLOw 

COMMON 



ENRG( 100) 

9 

RNG(IOO) 

COMMON 

FBOMBP (25,4) 

9 

FBOMBNf 25 ,4 ) 

9 

CPSP( 25,4) 

COMMON 

CPSN (25,4) 

9 

CNSP (25,4) 

9 

CNSN( 25,4) 

COMMON 

FNSP (25,4) 

9 

FNSN (25,4) 

9 

FBOMP (25,4) 

COMMON 

EBOMN ( 25 » 4 ) 

9 

EPROP (25,4 ) 

9 

EPRON (25,4) 

COMMON 

ENEUP (25,4) 

9 

ENEUN (25,4) 

9 

ENRGP (25) 

COMMON 

FNRGN ( 75 ) 

9 

XSMBP ( 25 ) 

9 

XSMBN ( 75 ) 

COMMON 

SNRG(IOO) 

9 

RBENROf 20 ) 

9 

Cl NRG (40) 

COMMON 

C2NRG ( 40 ) 

9 

CNRG ( 2 ) 

9 

SOFF( 100) 

COMMON 

RBF ( 20 ) 

9 

CnNKl(40) 

9 

C0NK2 ( 40 ) 

COMMON 

CONK ( 2 ) 

9 

LFNGTHf 8 ) 

9 

GMWT 

COMMON 

LSOFF 

9 

LRBF 

9 

LK 1 

COMMON 

LK2 

9 

LK 

9 

NOFCOM 

COMMON 

MOVE 

9 

KOSW{ 36 ) 

9 

KONSWT 

DIMFNSlON ENERGY! IDO) 

♦ D 

ATA(IOO) ,CONS 

T ( 2 ) 


EQUIVALENCE ( ENERGY (1) »ENRGP ( ] ) ) , (DATA (1 ) ,XSMBP(1 ) ) 

I ND»INDEX 

IF (IND .LT. 3) GO TD 50 
RATIO® 6.0231F-4/GMWT 
DO 100 J*l,2 

MX ■ 25* ( J-l ) +LFNG T H ( J+5 ) 

CONST! J)=DATA(MX )*RATIO 
KONSWT » KOSW( J+25 ) 

GO TO ( 100*108 1 »KONSW T 

108 MN* 25* ( J-l ) +1 
DO 109 K*MN,MX 

IF (FI^ERGY(K) .NE. 0.0) GO TO 30 
FNERGY(K) »1.F-10 
GO TO 31 

30 ENERGY ( K ) *AL0G1 0 ( FNFRGY ( K ) ) 

31 IF( DAT A ( K ) .NF. 0.0) GO TO 33 
DAT A ( K ) * l.F-10 

GO TO 109 

33 DATA(K)*ALOG10(DATA(K) ) 

109 CONTINUE 
100 CONTINUE 

GO TO 75 

50 FF=F 

MX « 2 5 * ( IND-1 )+LENGTH( IND+5) 

KONSWT * KOSW ( I ND+25 ) 

GO TO (52*51) .KONSWT 

51 EE= ALnGl 0 ( FF ) 

52 I F ( FF .LT. FNFRGYfMX ) ) GO TO 54 
XSFCT=cONST( IND) 

GO TO 75 

54 GO TO (56,55) ,IND 

55 IF(FF-FNRGN( 1 ) ) 110,120,56 

110 XSECT=0.0 
GO TO 75 

56 MN= 25*( IND-1 )+l 

CALL LAGRNGtFF, CROSS, FNFRGY(MN) ,DATA(MN) ,LFNGTH ( I ND+5 ) , 2 ) 

57 GO TO (59,5R) ,<ONSWT 



r> r\ n r> 


59 


58 • CR0SS=1 0.**CR0SS 
59 XSPTT=CRnSS*RATin 
75 RETURN 

120 CR0SS=XSMBN( 1 ) 

GO TO 57 
FND 

SIPFTC DOSFK LIFT ,DFCK 

SUBROUTINF DOSFX ( F ,VARB1 *VARB2 * I NOFX ) 

THIS SUBROUTINF CALCULATES (DE/DX ) =SOFF ,RBF ,NFUTRON FLUX TO DOSE 
CONVERSION FACTORS 
VARB1- DUMMY VARIARLF 
VARP2- DUMMY VARIABLE 


COMMON 

D2X ( 20 ) 

* 

MAX 

9 

FD ( 300 , 20 ) 

COMMON 

FI (300) 

♦ 

DPI (300 ) 

9 

PIRAR (300) 

COMMON 

OR ( 300 ) 

9 

OPPRM ( "*00 ) 

9 

N0D2X (20) 

COMMON 

X ( 20 ) 

9 

NOX 

9 

FNTOTS ( 200 ) 

COMMON 

DX ( 20 ) 

9 

PROPNOI 20 ) 

9 

Cl 

COMMON 

PDSBND 

9 

ndsbnd 

9 

BNDLOW 

COMMON 



ENRGt 100) 

9 

RNG( 100) 

COMMON 

EBOMBP (25,4) 

9 

EBOMBN ( 25 , A ) 

9 

CPSP ( 25 ,4 ) 

COMMON 

CPSN (25,4) 

9 

CNSP( 25,4) 

9 

CN?N( 25 ,4 ) 

COMMON 

FNSR ( 25 , A ) 

9 

FNSN (25,4) 

9 

FBOMP (25,4) 

COMMON 

FBOMN (25,4) 

9 

FPROR (25,4) 

9 

FPRON (25,4) 

COMMON 

FNEUP ( 25 » A ) 

9 

FNFUN (25,4) 

9 

RNRGP ( 25 ) 

COMMON 

ENRGN ( 75 ) 

9 

XSMBP ( 25 ) 

9 

XSMBN (75 ) 

COMMON 

SNRGtlOO) 

9 

RBFNRG (20) 

9 

C1NRG ( 40 ) 

COMMON 

C2NRGU0) 

9 

CNRG ( 2 ) 

9 

SOFE(IOO) 

COMMON 

RBF( 20) 

9 

C0NKK40) 

9 

C0NK2 (40) 

COMMON 

CONK (2) 

9 

LENGTH( 8 ) 

9 

GMWT 

COMMON 

LSOFE 

9 

LRBE 

9 

LK1 

COMMON 

LK2 

9 

LK 

9 

NOFCOM 

COMMON 

MOVE 

9 

K0SWI36 ) 

9 

XONSWT 


DIMENSION NOENTS ( 5 ) * ENERGY C 202 ) ,TABLE(202) *MAXVLU(5) ,VARB(5) , 
AEBND ( 2 ) *M I NS ( 5 ) 

DATA ( MINS { J ) * J = 1 *5)/l *101*121 *161 *201/ 

EQUIVALENCE (SNRG (1 ), ENERGY 1 1 ) ) , (SOFF ( 1 ) *T ABLE ( 1 ) ) * ( NOFNTS { 1 ) ♦ 
ALSOFE) , ( EBNDI 1 ) tPDSBND) 

REAL MAXVLM.NDSBND 

TF UNDEX.PT.il GO TO 10 
DO 10 J=1 * 5 

MX*MINS( J)+NOENTS( J)-l 
MAXVLU(J) « TABLE (MX ) 

KONSWT - KOSW ( J+31 ) 

GO TO (10,8) *KONSWT 

8 MIN*MINS(J) 

DO 9 K*M I N , MX 

IF ( ENERGY ( X 1 .NE. 0.0) GO TO 80 
ENERGY(K) « l.E-10 
GO TO 31 

30 FNRRGY ( K) *aL0G1 0 ( ENFRGY ( K1 1 

31 IF ( T ABLE ( X ) .NF. 0.0} GO TO 33 
TABLE(K)=1.E— 10 

GO TO 9 

83 TABLE (Kl*ALOG10(TAPLF(K) } 

9 CONTINUE 
10 CONTINUE 

GO TO 300 

100 IF (INDEX. GT. 3) GO TO 35 
INDM1 - INDEX-1 

150 IF (E.GT.EBND( INDM1 )) GO TO 154 
152 VARB1*0.0 
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* 


VARB2=0.0 
GO Tn 300 

154 MX = 2*1 NDM1 
MIN * MX -1 

404 DO 400 J»MIN,MX 
FE = F 

LIM=MINS( J ) +NOFNTS ( J ) — 1 
KONSWT * KOSW ( J+31 ) 

GO Tn (402,401 ), KONSWT 

401 EE * ALOGIO(FF) 

402 IF (FF-ENFRGY(LTM) ) 155.160,160 
160 VARB(J) = MAXVLUIJ) 

GO TO 400 

155 LIM=MINS(J) 

CALL LAGRNG (FF»VARB(J) .FNFRGY ( LIM ) ,TABLF(LIM) ,N0FNT5 ( J ) *2) 
GO Tn (400,403) .KONSWT 

405 VARB(J) = 10.**VARB(J) 

400 CONTINUE 

VARB1 » VARR(MIN) 

VARB2 * VARB ( MX ) 

500 RETURN 
550 MIN = 5 
MX = 5 
GO TO 404 
END 

SIPFTC RANGF LIST ,r>PCK 

SUBROUTINE RANGF (X.Y.INDFX) 

C RANGF-ENFRGY CALCULATIONS 


COMMON 

D2X ( 20 ) 

* 

MAX 

9 

FO ( 300,20 ) 

COMMON 

FI (500) 

» 

DFI (300) 

9 

FIBAR ( 300 ) 

COMMON 

OP (500) 

* 

OPPRM ( 5 00 ) 

9 

N0D2X ( 20 ) 

COMMON 

X ( 20 ) 

9 

NOX 

9 

FNTOTS ( 200 ) 

COMMON 

DX ( 20 ) 

9 

PR0PNn(20) 

9 

Cl 

COMMON 

PDSBND 

9 

NDSBND 

9 

BNDLOW 

COMMON 



FNRG ( 100) 

9 

RNG(IOO) 

COMMON 

FBOMBP (25,4) 

9 

FBOMBN (25,4) 

9 

CPSP( 25,4) 

COMMON 

CPSN (25,4) 

9 

CNSP (25,4) 

9 

CNSN (25,4) 

COMMON 

FNsP (25,4) 

9 

FNSN ( 25 ,4) 

9 

FBOMP (25,4) 

COMMON 

FBOMN (25,4) 

9 

FPROP (25,4) 

9 

FPRON (25,4) 

COMMON 

FNFUP (25,4) 

9 

FNFUN (25,4) 

9 

FNRGP(?^) 

COMMON 

FNRGN ( 75 ) 

9 

XSMSP ( 25 ) 


XSMBN (75) 

COMMON 

SNRG( ion) 

9 

RBENRG( 20 ) 

9 

C1NRG ( 40 ) 

COMMON 

C2NRG ( 40 ) 

9 

CNRG ( 2 ) 

9 

SOFE( 100) 

COMMON 

RBF ( 20 ) 

9 

C0NK1 (40) 

9 

CONK 2 (40) 

COMMON 

CONK ( 2 ) 

9 

LENGTH ( 8 ) 

9 

GMWT 

COMMON 

LSOFE 

9 

LRBE 

9 

LK1 

COMMON 

LK2 

9 

LK 

9 

NOFCOM 

COMMON 

MOVF 

9 

K0SWI36 ) 

9 

KONSWT 

FQUTVALENCF (KnSW<25) 

.K0SW25 ) 



DIMENSION CON ( 2 ) 





IF (INDEX ,GT » 1) Gn 

TO 

20 




CON ( 1 ) *ENRG ( 1 ) 

CON ( 2 ) *RNG ( 1 ) 

Ll-LENGTH(l) 

GO TO (5o,ir) ,X0SW25 

10 DO 1? 1*1, LI 

IF ( FNRG ( I ) , NF« 0 , 0 ) on Tn 11 
FNRG(I) « l.F-10 
GO Tn 12 

11 FNRG(I)*AL0G10(FNRO(D) 

12 IF (RNG(I) , NF, 0.0) GO TO 14 


0010 

0020 
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RNG(I) * l.F-lo 
60 TO 13 

14 RNG(I)= AL0610 (RNG (T1 I 
13 CONTINUE 
GO TO 50 
20 XX = X 

IMl=INDEx-l 

I F ( XX .GF. CON < IM1 ) ) GO TO 25 

y*o,o 

GO TO 50 

25 GO TO (27*26) ,K0SW25 

26 XX=AL0G10(XX) 

27 GO TO (25,29) *IM1 

25 CALL LAGRNG I XX ,Y ,FNRG ,RNG,L1 ,2 ) 
GO TO 30 

29 CALL LAGRNG ( XX,Y,RNG,FNRG,Ll ,2 ) 

30 GO TO ( 50,31) ,K0SW25 

31 Y=10.*»Y 
50 RETURN 

END 

$ IBFTC SORT LIST, DECK 

SUBROUTINE SORT ( E ,FRNO ,MAXX, I ) 


DIMENSION ERNG ( 300 ) 

DO 10 J=1,MAXX 0060 

IF ( FRNG ( J ) -F ) 14,12,10 0070 

10 CONTINUE 0050 

12 I=J 0090 

GO TO 15 0100 

14 I=J-1 0110 

15 RETURN 0120 

END 0130 

$ IBFTC YIFLDS LIST ,nFCK 


SUBROUT INF YIELDS ( , ANS,NOC, INDFX ) 

C CALCULATION OF YIELDS AS A FUNCTION OF ENERGY OF BOMBARDING 

C PARTICLE. 

C 


COMMON 

D2X ( 20 ) 

9 

MAX 

9 

FO ( 300 , 20 ) 

COMMON 

El (300) 

9 

DEI (300) 

9 

EIBAR (300) 

COMMON 

OP ( 300 ) 

9 

OPPRM ( 300 ) 

9 

N0D2X ( 20 ) 

COMMON 

X ( 20 ) 

9 

NOX 

9 

ENTOTS ( 200 

COMMON 

DX ( 20 ) 

9 

PROPNO ( 20 ) 

9 

Cl 

COMMON 

PDsBND 

9 

NDSBND 

9 

BNDLOW 

COMMON 



FNRG(IOO) 

9 

RNG(IOO) 

COMMON 

EBOMBP (25,4) 

9 

FB0M8N( 25,4) 

9 

CPSP( 25,4) 

COMMON 

CPSN (25,4) 

9 

CNSP ( 25 ,4) 

9 

CNSN (25,4) 

COMMON 

ENSP (25,4) 

9 

ENSNt 25,4) 

9 

EBOMP (25,4 

COMMON 

FBOMN (25,4) 

9 

EPROP (25,4) 

9 

EPRON (25,4 

COMMON 

ENEUP (25,4) 

'* 

FNFUN (25,4) 

9 

ENRGP ( 25 ) 

COMMON 

ENRGN ( 75 ) 

• 

XSMBP ( 25 ) 

9 

XSMBN ( 75 ) 

COMMON 

SNRG(IOO) 

♦ 

RBENRG ( 20 ) 

9 

C1NRG (40) 

COMMON 

C2NRG ( 40 ) 

9 

CNRG( 2 ) 

9 

SOFF (100) 

COMMON 

RBE ( 20 ) 

t 

CONK! (40) 

9 

C0NK2 (40) 

COMMON 

CONK ( 2 ) 

* 

LENGTH! 5 ) 

9 

GMWT 

COMMON 

LSOFE 

• 

LRBE 

9 

LK1 

COMMON 

LK2 

9 

LK 

9 

NOFCOM 

COMMON 

MOVE 

* 

KOSW( 36 ) 

9 

KONSWT 


DIMENSION ANS(3) ,EB0MB(20O) ,CPS(200) ,CNS(200) ,FNS(200) 

EQUIVALENCE (EBOMB(l) ,FB0MBP(1,1 ) ) » ( CPS ( 1 ) ,CPSP( 1 ,1 ) ) , 

A ( CNS ( 1 ) ,CNSP ( 1 , 1 ) ) , ( FNS ( 1 ) ,ENSP (1,1)) 

CPS ( I ) -CASCADE PROTONS, CNSU J-CASCADE NFUTRONS, ENS ( I ) -EVAPORATION 

C NEUTRONS 
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IF (INDEX ,GT. 1 ) GO TO 100 

00 321 <*1,2 

ICON SWT * KOSWOC+27) 

GO TO (321,40) ,KONSVfT 
60 DO 320 L=1 .NOFCOM 

MN=100*(K-1)+25*(L-1)+1 
MX=MN+LENGTH ( K+l ) -1 
DO 320 I =MN . MX 

IF (FBPMB(I) .NF. 0.0) GO TO 1 
EBOMB(I) * l.E-10 
GO TO 2 

1 EBOMB(I) * ALOG10(FB0MR(l) ) 

2 IF CCPS(I) .NF. 0 . 0 ) r,n TO 3 
CPS(I) * l.E-10 

GO TO 6 

3 CPS ( I ) *AU0G10 ( C°S ( I ) ) 

6 IF (CNS(I) .NE. 0.0) GO TO 5 
CNS(I) = l.E-10 
GO TO 6 

*5 CNS ( I ) * AL0G10 ( GNS ( I ) ) 

6 IF (FNS(I) .NF. 0 . 0 ) op TP 7 
ENS ( I ) = l.F-in 

GO TP 320 

7 ENS ( I ) = ALOG1 0 ( FNS ( I ) ) 

320 CONTINUE 

321 CONTINUE 
GO Tp 200 

100 F=FF 

XONSWT = X0FW( INDEX+26 ) 

GO TO ( 106,102 ) .KONSWT 
102 F=AL0G10(F) 

104 MN*100*( INDFX-2)+25*(N0C-l ) + l 

CALL LAGRNG(E.ANSd) .EBOMB(MN) ,CPS(MN) .LENGTH ( INDFX ) .2) 

CALL LAGRNG ( F » ANS ( 2 ) .EBOMB(MN) ,CNS(MN) .LfNGTH t INDEX ) ,2) 

CALL LAGRNG(F,ANS(3! .FPOMB(MN) ,FNS(MN) .LFNGTH ( INDEX ) ,2) 

GO TO ( 200, iofl ) .KONSWT 
108 DO 110 J*1 , 3 
110 ANS(J)»10.**ANS( J) 

200 RETURN 
END 

SIBFTC CASNRG LIST ,PFC< 

SUBROUTINE CASNRG ( FF. ANS ,NOC » TNOFX ) 

CALCULATES PNFRGY OF CASCADE PROTONS AND NEUTRONS AS A FUNCTION 

C .... .OF THF FNFRGY OF THF BOMBARDING PART I CLF . 


COMMnN 

D2X( 20) 

* 

MAX 

9 

FT ( 300,20 ) 

COMMON 

FI (300) 

9 

DFT (30P) 

9 

•dBAR ( 300 ) 

COMMON 

OP { 300 ) 

9 

OPPRM ( 300 ) 

9 

NGD2X ( 20 ) 

COMMON 

X ( 20 ) 

9 

NOX 

9 

BNTOTS ( 200 ) 

COMMON 

DX ( 2 0 ) 

9 

PROPNO ( 20 ) 

9 

Cl 

COMMON 

PDSBND 

# 

NDSBND 

9 

BNDLOW 

COMMON 



FNRG(IOO) 

9 

RNG( 100 ) 

COMMON 

EBOMBP (25.6) 

9 

FBOMBNI 25 .6) 

9 

CPSP( 25,6) 

COMMON 

CPSNI25.6) 

9 

CNSP (25,6) 

9 

CNSN( 25,6) 

COMMON 

ENSP (25,6) 

9 

FNSN (25.6) 

9 

FBOMP (25,4) 

COMMON 

EBOMN (25.6) 

9 

EPROP (25,6) 

9 

EPRON (25,4) 

COMMON 

ENEUP (25.6) 

9 

FNFUN (25,6) 

9 

FNRGP ( 25 ) 

COMMON 

FNRGN ( 75 ) 

9 

XSMBP ( 25 ) 

9 

XSMBN ( 75 ) 

COMMON 

SNRG(IOO) 

9 

RBFNROI 20) 

9 

C1NRG ( 40 ) 

COMMON 

C2NRG ( 40 ) 

9 

CNRG ( 2 ) 

9 

SOFF( 100) 

COMMON 

RBF( 20) 

9 

C0NK1 (40) 

9 

C0NK2 (40) 

COMMON 

CONK ( 2 ) 

9 

LFNGTH! 8 ) 

9 

GMWT 
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COMMON 

LSOFE 

, LRBE 

, LK1 

COMMON 

LK2 

, LK 

, NOFCOM 

COMMON 

MOVE 

, KOSW( 36 ) 

, KONSWT 

DIMENSION ANS ( 2 ) 

, FBOMB (200), FPRO ( 200) 

, FNFU ( 200 ) 


EQUIVALENCE (EB0MB(1 ) ,FBOMP( 1 ,1 ) 1 , (EPROd ) ,FPROP ( ) ,1 ) ) , 

A (FNFU(l).PNFUP(l.in 

C FBOMB- ENERGY nF BOMBARD T NG PAPTICLF 

C FPRO- SECONDARY PROTON ENERGY 

C FNFU- SECONDARY NEUTRON ENERGY 

IF (INDEX .GT. II GO TO 200 
DO 121 K=l,2 
KONSWT = K0SW(K+29) 

GO TO ( 121,110) .KONSWT 
110 DO 120 L*1,N0FC0M 

MN=lnn*(K-l)+25*(L-l 1+1 
MX=MN+LENGTH ( K+3 1 -1 
DO 120 J=MN,MX 

IF ( EB0MB( J) .NF. 0.0) on TO 1 
EBOMB(J) = l.E-10 
GO TO 2 

1 FBOMB ( J)=AL0G10f FBOMB( J) I 

2 I F ( FPRO ( J I .NF. 0.0) GO TO 1 
FPRO(J) = l.F-10 

GO TO 4 

3 EPRO ( J ) “AL0G1 0 ( FPRO ( J ) ) 

4 IF(ENEUIJ) .NF. 0.0) GO TO 5 
ENFU(J) * l.E-10 

GO TO 120 

5 FNEU( J)=AL0G10(FNFU( J)) 

120 CONTINUE 

121 CONTINUE 
GO TO 300 

200 E=EF 

KONSWT = KOSW( INDEX+28 ) 

GO TO (204,202) .KONSWT 
202 E-ALOGIO(F) 

204 MN*1 00* ( I NDFX-2 ) +28* ( NOC— 1 )+l 

CALL LAGRNG(F,ANS(1 ) .FBOMB ( MN ), FPRO ( MN ) .LFNGTH ( I NDEX+2 ) . 2 ) 

CALL LAGRNG(F,ANS(2) .FROMB(MN) .FNFU(MN) .LFNGTH ( TNDFX+2 ) ,2) 

GO TO (300,206) .KONSWT 
206 ANS ( 1 ) =10«0**ANS ( 1 ) 

ANS ( 2 ) *10.0** ANS ( 2 ) 

300 RETURN 
END 

SIBFTC LAGRNG LIST.DFCK 

SUBROUTINE LAGRNG ( XX , YY ,XTAB ,YT AB »L IMT T .RANK ) 

INTERPOLATION SUBROUTINE BASFD ON LAGRANGE-S FUNDAMENTAL FORMULA 

C FOR INTERPOLATION 

DIMENSION XTAB(LIMIT) .YTABCLIMIT) ,DIFS(15) 

INTEGER ORDER, HALF, 0RDM1, RANK 

x*xx 

MAXNO=LIMIT 

ORDER=RANK 

DO 10 INDFX*10,MAXN0,1P 
IF (INDEX .GF. MAXNO) GO TO 15 
IF (X - XTAB ( I NDFX ) ) 15,55,10 

10 continue 
15 J= I NDEX-9 

DO 25 INDEX*J, MAXNO 
IF , ( X-XTAB ( I NDFX ) ) 28,55,25 
25 CONTINUE 
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28 HALE=0PDER/2 
0RDM1*0R0ER-1 

IF (INDEX .GT. HALF+1I GO TO 33 

MIN*1 

GO TO 39 

33 IF (INDEX .LT. MAXNO-HALF+1 ) GO TO 37 
MIN*MAxNO-ORDM1 
GO TO 39 

37 MIN*INDEX-HALE 
39 MAx=MI N+ORDM 1 
MINMl*MlN-l 

g DO 42 J»l, ORDER 

O I NDEX*MI NM1+J 

*? 42 DIFS ( J ) = X— XTAB ( I NDFX ) 

w Y = 0.0 

DO 50 J*MIN,MAX 
TFRM* YTAB(J) 

DO 45 INDFX = MIN»MAX 
IF (J .EQ. INDFX) GO TO 45 

44 MARX* INDFX-MINM1 

TERM = TFRM*DIFS(MARK)/(XTAB(J)-XTAB( INDFX) ) 

45 CONTINUE 
50 Y*Y+TFRM 

yy*y 

52 RETURN 
55 YY*YTAB( INDFX) 

GO TO 52 
END 

SIBFTC PROPTY LIST »DFCK 

SUBROUTINE PROPTY ( INDFX .LAYER ) 

THIS subroutine TRANSFFRS the material property data from tape 3 
TO TAPE 4 ( OR DISC STORAGE) FOR LATFR "SE. THE TABLES OF FLUX TO 
DOSE CONVERSION FACTORS ARE TRANSMITTED FROM TA°F 3 TO CORF 
STORAGE. AT THF APPROPRIATE TIME THE PROPERTY DATA FOR THE CHOSEN 
MATERIAL IS TRANSFERRED FROM TAPE 4 ( OR DISC STORAGE) TO CORE AND 
ANY SUBROUTINES USING THIS PROPERTY DATA ARF INITIALLED. 


COMMON 

D2X ( 20 ) 

* 

MAX 

9 

F0< 300,20) 

COMMON 

El (300) 

9 

DEI (300) 

9 

EIBAR (300) 

COMMON 

OP (300 ) 

♦ 

OPPRM (300 ) 

9 

N 0 D 2 X ( 20 ) 

COMMON 

X ( 20 ) 

9 

NOX 

9 

ENTOTS ( 200 ) 

COMMON 

DX ( 20 ? 

9 

PROPNO ( 20 ) 

9 

Cl 

COMMON 

PDSBND 

9 

NDSBND 

9 

BNDLOW 

COMMON 



ENRGt 100) 

9 

RNG(IOO) 

COMMON 

EBOMBP (25.4) 

9 

EBOMBN( 25*4) 

9 

CPSP( 25,4) 

COMMON 

CPSN (25*4) 

9 

CNSP( 25*4) 

9 

CNSN( 25,4) 

COMMON 

ENSP(?5.4) 

9 

ENSN (25*4) 

9 

EBOMP (25,4) 

COMMON 

EBOMN ( 25 »4 ) 

9 

EPROP (25.4) 

9 

EPRON (25,4) 

COMMON 

FNFUP (25*4) 

9 

FNEUNI 25.4) 

9 

ENRGP ( 25 ) 

COMMON 

FNRgN ( 75 ) 

9 

XSMBP ( 25 ) 

9 

XSMBN ( 75 ) 

COMMON 

SNRG(IOO) 

9 

RBENRGt 20) 

9 

C 1 NR G ( 4 0 ) 

COMMON 

C2NRG ( 40 ) 

9 

CNRGI2) 

9 

SOFE(IOO) 

COMMON 

RBF ( 20 ) 

9 

CONKl (40) 

9 

C0NK2 (40) 

COMMON 

CONK ( 2 ) 

9 

LENGTH! 8 ) 

9 

GMWT 

COMMON 

LSOFE 

9 

LRBE 

9 

LX 1 

COMMON 

LX2 

9 

LX 

9 

NOFCOM 

COMMON 

MOVP 

9 

K0SWO6 ) 

9 

XONSWT 


FQUIVALENCF ( XOSVM 1 3 ) .K0SW1 3 ) 

EQUIVALENCE ( LFNGTH ( 1 ) ,L1 ) . (LFNGTH(2) .L2 ) , ( LFNGTH ( 3 1 *L3) * (LFNGTH 
A ( 4 ) ,L4 ) » ( LFNGTH ( 5 J »L5 ) , ( LENGTH ( 6 ) »L 6 ) * ( LENGTH ( 7 ) *L7 ) . ( LFNGTH ( 8 ) » 
BL 8 ) 

EQUIVALENCE (LFNTH(l) ,LX1) »(FTDC0N(1) .CONXl(l) ) 
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■DIMENSION LFNTH ( 3 ) ,FTDCON(40.25 
INTEGER PROPNO 
L=L AYER 

IF (INDEX .GT. 1) GO TO 100 
REWIND 4 
LASTN0=0 
I SAVE= 1 

READ (3) NOMATl 
DO 10 N=1»NDMAT1 

READ (3) NO, NOFCOM, GMWT.L1,L2,L3,L4,L5.L6,L7, (ENRG( J) ,RNG( J) ,J=1. 
ALl ) , < (FBOMP( J.K) ,EPROP( J.K) ,FNFUP( J.K) ,J=1 ,L2) ,K = 1 .NOFCOM) , ( (EBOMN 
B( J.K) ,EPRON( J.K) ,ENEUN( J.K) , J=l ,L3 ) ,K=1, NOFCOM) , ( ( FBOMBP ( J , K ), CPSP 
CC J.K) .CNSP( J.K) *ENSP (J.K) .J=1»L4> .K=l .NOFCOM) , ( (EBOMBN( J.K) ,CPSN( 
DJ.K) ,CNSN( J.K) .ENSN( J.K) ,J*1 ,L5) »K=1, NOFCOM) . ( ENRGP ( J ) » XSMBP ( J ) »J= 
El »L6 ) . ( ENRGN ( J ) . XSMBN ( J ) , J=1 ,LT ) 

10 WRITE (4) NO, NOFCOM, GMWT »L1 » L 2 »L3 »L4 »L5 .L6.L7, ( FNRG ( J ) ,RNG ( J ) ,J=1. 
ALl > , ( (EBOMP( J.K) , EPROP (J.K) ,FNFUP( J.K) ,J=1 ,L2 ) ,K=1 , NOFCOM) , ( ( EBOMN 
B( J.K) ,fPRON( J.K) ,ENFUN( J.K) , J= 1,L3 ) ,K*1 .NOFCOM) , ( ( FBOMBP ( J ,K ) ,CPSP 
C( J»K) ,CNSP( J.K) ,ENSP( J.K) ,J=1 ,L4) ,K=1 .NOFCOM) , ( ( EBOMBN (J.K) .C^SNf 
DJ.K) .CNSN(J.K) ,FNSN( J.K) , J=1 , L5 ), K = l , NOFCOM ),( FNRC,p ( J ), XSMBP ( J ) ,J» 
El ,L6) , (ENRGN ( J) .XSMBN ( J ) ,J*1,L7) 

REWIND 4 

READ (3) LRBE»LK1 .LK2.LK, (RBENRr,( J) ,RBE( J) ,J*1 ,LRBF) , ( Cl NRG ( J ) , 
1C0NKK J) ,J=l,LKl ) , ( C?NRG ( J ) ,C0NK2( J) ,J = 1,LK?) , { CNRG ( J ) .CONK ( J ) , 
2J*1.LK) 

IF (K0SW13 .FO. 2) GO TO 200 

DO 210 J*1 ,2 

LIMT=LFNTH(J) 

DO 210 K*1,LIMT 

210 FTDC0N(K,J)*FTDC0N(K,J)*360O, 

200 READ (3) NOMAT 2 
DO 12 N“1 .N0MAT2 

READ (3) N0,LS0FE,(SNRG(J),S0FE(J) .J-l.LSOFF) 

IF (NO .EG. L ) GO TO 20 

12 CONTINUE ? 

WRITE (6,14) L 

14 FORMAT! 1H0RX , 5SHDATA TAPE DOES NOT CONTAIN DE/DX TABLE FOR MATERIA 
1L NUMBERI4) 

16 REWIND 3 
STOP 

20 REWIND 3 

CALL DOSFK ( DUMMY .DUMMY .DUMMY, 1 ) 

RETURN 

100 !F(PROPNO(L) .GT. LASTNO ) GO TO 105 
MAVG« ( PROPNO ( L )+LASTNO ) /2 
IF( MAVG .GT. 100) GO TO 110 
REWIND 4 
I 5AVE=1 
GO TO 105 

110 KOUNT«LASTNO-PROPNO(L) 
isavf«isavf-kount 
KOUNT «KOl.)NT+l 
DO 103 LL*1, KOUNT 
103 BACKSPACE 4 
105 DO 120 NN* I SAVE. NOMAT 1 

READ (4) NO, NOfcOM, GMWT,L1,L2, L3, L4, L5, Lfe, L7, ( ENRG( J ) ,RNG(J) ,J*1, 
ALl ) , ( (FBOMP( J.K) .EPROP (J.K) ,ENFUP( J.K) ,J=1 ,L2 ) »K*1 .NOFCOM) , ( (EBOMN 
B( J.K) ,EPRON( J.K) ,ENFUN( J.K) , J»1,L3 ) ,K*i .NOFCOM ),(( FBOMBP ( J , K ) ,CPSP 
C( J.K) ,CNSP( J.K) ,ENSP( J.K) ,J*1,L4) ,K*l .NOFCOM) , ( ( FBOMBN ( J.K) ,CPSN( 
DJ.K) ,CNSN( J.K) ,FNSN( J.K) ,J*1 »L5 ) , K*1 .NOFCOM ) , ( ENRGP ( J ) , XSMBP ( J ) ,J= 
El »L6) , ( ENRGN ( J) .XSMBN ( J) ,J=1»L7) 
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IF ( PROPNO (L) .EQ. NOl 60 TO 150 
120 CONTINUE 

WRITE (6*4) PROPNO(L) 

4 FORMAT ( 1H08X .51HPR0GRAM CANNOT FIND DATA TABLES FOR MATERIAL NUMBE 
1RI4) 

STOP 

150 LASTNO =PROPNO(L) 

ISAVE-NN 

CALL RANGE (DUMMY,DUMMY*1) 

IF (INDEX .LT. 3) GO TO 151 
IF ( K0SW13 .EQ. 1) C1=C1*3600. 

CALL XS (DUMMY, DUMMY *3) 

CALL YIELDS ( DUMMY, DUMMY, NOFCOM, 1 ) 

CALL CASNR6 ( DUMMY .DUMMY .NOFCOM , 1 ) 

151 RETURN 
END 

SIBFTC FLUXFQ LIST ,PFCK 


SUBROUTINE FLUXEQ (E,FLUX,NO) 

C CALCULATES INITIAL INCIDENT PROTON SPECTRUM AS A FUNCTION OF 

C INITIAL INCIDENT ENFROIFS. 

C 

C IF NO FOUALS- flux EQUAL S- 

C 1 A*F**(-B), 

C 2 A*EXP ( — P ( F) /PP) ( N ( GREATER THAN P)), 

C 3 TABLF OF FLUX VS. E AND INTERPOLATION, 

C 4 A ( F) *EXP ( — B ( E ) ) , 

C 5 10.0**( A1+A2* C +A3*F**2+A4*F**3) , 

C 6 10.0** (A1+A2*L0G(E) +A3* ( LOG ( F ) ) **2+A4* ( LOG ( F) ) **3 ) , 

C 7 -A/P0*FXP(-C1»P(E)/P0)*P1(F)/D(F) 

C WHERE E IS THE GIVEN INCIDFNT FNERGY. 


COMMON 

D2X ( 20 ) 

♦ 

MAX 

9 

EO ( 300,20 ) 

COMMON 

El (300) 

♦ 

DFI (300) 

9 

FIBAR (300) 

COMMON 

OP (300 ) 

9 

OPPRM (300 ) 

9 

N0D2X ( 20 ) 

COMMON 

X ( 20 ) 

9 

NOX 

9 

ENTOTS ( 200 ) 

COMMON 

DX ( 20 ) 

9 

PROPNO (20) 

9 

Cl 

COMMON 

PDSBND 

9 

NDSBND 

9 

BNDLOW 

COMMON 



FNRG( 100) 

9 

RNG( 100 ) 

COMMON 

EBOMBP (25,4) 

9 

EBOMBN ( 25,4) 

9 

CP SP( 25,4) 

COMMON 

CPSN( 25,4) 

9 

CNSP( 25,4) 

9 

CNSN( 25,4) 

COMMON 

ENSP (25,4) 

9 

ENSN( 25,4) 

9 

EBOMP (25,4) 

COMMON 

EBOMN (25,4) 

9 

EPROP (25,4) 

9 

EPRON (25,4) 

COMMON 

FNFUP< 25,4) 

9 

FNFUN (25,4) 

9 

FNRGP (25) 

COMMON 

ENRGN ( 75 ) 

9 

XSMBP ( 25 ) 

9 

XSMBN ( 75 ) 

COMMON 

SNRG(IOO) 

9 

RBFNRG ( 20 ) 

9 

C1NR6 ( 40 ) 

COMMON 

C2NRG ( 40 ) 

9 

CNRG ( 2 ) 

9 

SOFF( 100) 

COMMON 

RBE( 20 ) 

9 

C0NK1 (40) 

9 

CONK 2 (40) 

COMMON 

CONK ( 2 ) 

9 

LENGTHt 8 ) 

9 

GMWT 

COMMON 

LSOFE 

9 

LRBE 

9 

LK1 

COMMON 

LK2 

9 

LK 

9 

NOFCOM 

COMMON 

MOVE 

9 

KOSW( 36) 

9 

KONSWT 


DIMENSION FFFF(IOO) ,PR0TS(100) ,A(4) ,B(4) 

EE=E 

IF (M0VE.E0.2 ) GO TO ( 10,20*30,40,50,60,70) ,N0 
MOVE * 2 

GO TO (1, 2, 3, 4, 5, 6, 7) , NO 
1 READ (5*100) A ( 1 ) »B ( 1 ) 

100 FORMAT ( 4E1 2 • 5 ) 

WRITE (6,102) A ( 1 ) ,B ( 1 ) 

102 FORMAT ( lH08X,22H0 = A*F**(-B) WITH A =1PF13.5,8H AND B «F13.5/1H+ 
A8X.1HI ) 
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10 'PHI * All )*FE**(-B(1 ) ) 

1000 FLUX=PH I 
2000 RETURN 

2 READ <5*100) A<1 ) *P0 
WRITE (6,202) A ( 1 ) *P0 

202 F0RMAT(1H0RX,44H0(GREATER THAN P) = A*FXP(-P( F ) /P0) WITH A =1PF13. 
A5,9H AND PO =F1 3 » 5/1 H+8X ,1HI ) 

20 P = 938. 26*SQRT( (EF/938. 26+1. 0)**2-1.0j 
PHI * A<1 )*FXP(-P/PO) 

GO TO 100 

3 WRITE (6,302) 

302 FORMAT ( 1H08X43H0 CALCULATED FROM TABLF Of FLUX VS. ENFRGY. / 1H+8X , 
A1HI ) 

RFAD (5*303) NOfNTS. (FFRF( I ) ,PROTS( I ) * 1 = 1 .NOFNTS) 

303 FORMAT ( IA/(FR.0,E10.3,F8.0,E10.3*F8.0,E10.3,FR.o,E10.3 ) ) 

DO 306 1 = 1 .NOENTS 

EEEF( I )=AL0G10(FEEF( I ) ) 

306 PROTS ( I ) =AL0G1 0 ( PROT S ( I ) ) 

30 EE=AL0Gl0 ( EE ) 

CALL LAGRNGf EE ,PHI ,EEEE , PROTS ,NOENTS , 2 ) 

PH 1=10, **PHI 
GO TO lOn 

A READ (5*100) A ,B 

WRITE (6*402) ( I » A( I ) , I =1 ,4 ) , ( I ,B ( ! ) , !=1 ,4 ) 

402 FORMAT ( 1H08X*91H0 = A ( F ) *E** ( -B ( E ) ) WITH A ( E ) AND B ( F ) OF THF FORM 

1 C ( E ) = Cl + C2*E + C3*E**2 + C4*E**3 AND/1H+8X , 1HI /9X , 4 ( 4X » 1HAI 1 , 

2 2H = 1PE13.5.1H* ) »4H AND/9X , 4( 4X , 1HB 1 1 * 2H = IDE 13 . 5 » 1H . ) ) 

40 SUMB * B ( 4 ) 

DO 42 1*1*3 
11=4-1 

42 SUMB = 5UMB*EE+B (IT) 

50 PHI = A ( 4 ) 

DO 46 1*1*3 
II = 4-1 

46 PHI = PHI*EE+A(II) 

IF (NO - 5) 48,5060*5060 
48 PHI =PHI *EXP t -SUMB ) 

GO TO 100 
5060 PHI =10, **PH! 

GO TO 100 

5 READ (5*100) A 

WRITE (6,502) ( I ,A( T ) ,1*1,4) 

502 FORMAT (1H08X,42HL0G 0 = A1 + A2*E + A3*E**2 + A4*E**3 WITH/1H+12X, 
A1HI /13X»4(4X,lHAI 1 ,2H *1PE13 . 5 ,1H , ) ) 

GO TO 50 

6 READ (5,100) A 

WRITE (6,602) d , A ( I ) , 1*1 ,4 ) 

602 F0RMAT(1H08X,61HL0G 0 = A1 + A2*L0G(E) + A3* ( LOG ( F ) ) **2 + A4*(L0G( 
1E))**3 WI TH/1H+1 2X,1HI/13X,4(4X, 1HAI 1 , 2H =1PE 1 3 . 5 , 1H , ) ) 

60 EE * AL0G1O ( EE ) 

GO TO 50 

7 READ (5*100) A(1 ) ,P0 
WRITE (6,702) A ( 1 ) ,D0 

702 F0RMAT(lH0BX*53H0(-DN/DF) = A /PO*FXP ( -C1*P ( F ) /Do ) *Pl ( F ) /P ( F ) WITH 
1A =1PF13.5,9H AND PO =D13.5/1H+8X *1HI ) 

70 PI = FF/938.26 + 1.0 
P = SORT ( P1*P1-1 • 0 ) 

PHI = A ( 1 )/P0*FXP(-93B.26*P/P0)*Pl/P 

GO TO 100 

END. 

SIBFTC INVALU LIST, DECK 
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SUBRDUTINF INVALU 




INITIAL 

DATA CALCULATIONS 



MOVF = 

1, INTEGRAL SPECTRUM 



MOVE = 

2, DIFFERENTIAL 

SPECTRUM 



COMMON 

D2X( 20) 

, max 

* 

FO ( 300,20) 

COMMON 

El (300) 

• DEI (3005 

* 

FIBAR (300) 

COMMON 

OP ( 300 ) 

, OPPRM (300) 

9 

N0D2X ( 20 ) 

COMMON 

X ( 20 ) 

, NOX 

9 

FNTOTS ( 200 ) 

COMMON 

DX( 20) 

, PROPNO ( 20 5 

9 

Cl 

COMMON 

PDSBND 

, NDSBND 

9 

BNDLOW 

COMMON 


ENRG( 100) 

9 

RNG(IOO) 

COMMON 

EBOMBP (25*4) 

, EBOMBN (25,4) 

9 

CPSP ( 25,4) 

COMMON 

CPSN (25*4) 

. CNSP (25,4) 

9 

CNSN( 25,4) 

COMMON 

ENSP (25*4) 

, ENSN (25,45 

9 

EBOMP (25,4) 

COMMON 

EBOMN (25*4) 

, FPROP ( 25,4) 

9 

FPRON (25,4) 

COMMON 

FNFUP (25,4) 

* FNFUN (25,4) 

9 

FNRGP (25) 

COMMON 

FNRGN ( 75 ) 

, XSMBP ( 25 ) 

9 

XSMBN (75) 

COMMON 

SNRG( 100) 

, RBFNRG ( 20) 

9 

C1NRG ( 40 ) 

COMMON 

C2NRG ( 40 ) 

, CNRG ( 2 ) 

9 

SOFE( 100) 

COMMON 

RBE ( 20 ) 

, C0NKK40! 

9 

C0NK2 ( 40 ) 

COMMON 

CONK ( 2 ) 

, LENGTH ( 8 ) 

9 

GMWT 

COMMON 

LSOFE 

, LRBE 

9 

LK 1 

COMMON 

LK2 

, LK 

9 

NOFCOM 

COMMON 

MOVIE 

, KOSW( 36 ) 

9 

KONSWT 

EQUIVALENCE (KOSW( 3), 

K0SW3 ) , ( KOSW( 1 

1 ) 

,K0SW1 3 ) * ( KOSW( 15 ) ,K0SW15) 

A 

( KOSW( 19 ) , 

K0SW19 ) 



DIMENSION NO I NTS (2) *E0MAX(25.2) ,N0 1 NCR ( 25 , 2 1 *KNTR(2) ,ESPEC ( 200 * 2 ) 


ADELTF(25,2),K0UNT(20),TITLE(m ,01 ST (300) 

EQUIVALENCE (EIBARtl ) ,DIST( 1 ) ) 

integer propno.eqno 
D xm*xm 

D2X(1)=DX(1) /FLOAT ( N0D2X ( 1 ) ) 

IF(NOX-l) 5,9,5 
5 DO 7 J=2 ,NOX 
C DX-LARGE DELTA X 

C X(J)- PRINT BOUNDARY MEASURED NORMAL TO INCIDENT FACE 

DX< J)=X(J)-X( J-l ) 

C D2X ( J ) -SMALL DFLTA X 

c N 0 D 2 X ( J ) -Number of small delta x in large delta x 

7 D2X(J)=DX(J)/FL0AT(N0D2X(J) ) 

9 NOX D l=NOx+l 

GO TO (250,260) ,K0SW19 
260 READ (5*262) KEI * (FI ( J) , J*1 *KFI) 

262 FORMAT( I3/(8F9.0) ) 

GO TO A 5 
250 KLIM-1 

21 READ (5,22) NOI NTS ( KLIM ) 

22 FORMAT (13) 

INTNO=NOlNTS(KLTM) 

READ (5*24) ( EOMAX ( L *KL I M) , NO J NCR (L *KL IM ) ,L=1 , I NTNO ) 

24 FORMAT ( F8 ,0 . 1 4 , F8 . 0 . 14 *F8 . 0 , 14 , F8 .0 , 1 4 ,F8 .0 * 14 ,F8. 0 . 1 4 ) 
IF (KLIM .EQ. 2) CO Tn 228 
220 GO TO (225,222) .K0SW15 
222 KLI M=2 
GO TO 21 

225 READ (5*9999) NDF 
9999 FORMAT ( 1 3 ) 

ENDF=NDE 

228 DO 2011 L = 1 * KL I M 
EL0WER=0« 

KONTUR *1 


0010 

0020 

0030 

0040 


0160 

0170 

0180 

0190 


0200 


0210 

0220 
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•ESPEC(l»L)=n.O 

INTN0=N0INT?(L) 

DO 201 J=1,INTN0 

r DELTF ( J»L ) -SMALL DFLT A F AT Fx!T FAFF 

DELTF ( J,L)=(EOMAX( J,L)-ELOWER ) /FLOAT ( NOI NCR (J ,L ) ) 

ELOWFR=EOMAX( J,L) 

LIMIT=NOl NCR ( J *L ) 

DO 1202 1*1, LIMIT 
KONTUR = KONT(jR+l 

1202 FSPFC ( KONTUR ,L ) =ESPFC (XONTUR— 1 ,L )+DFLTF ( J,L ) 

201 ESPFC(XONTUR,L)=FOMAX(J,L) 

2011 XNTP(L)=XONUlR 

XLI V=1 

DO 10 J = 1 »NOX 
I=N0XP1-J 

IF (J . NE. 1) 00 TO 12 

15 READ (5,16) FIMAX 

16 FORMAT (F12.5) 

FI ( 1 ) = FIMAX 

DO 400 11=1, NOX 

IF( II .EO. 1 ) GO TO 405 

402 IF (PROPNO(II) .EQ, PROPNO ( I I -1 ) ) GO TO 404 

403 CALL PROPTY ( 2 , 1 1 ) 

CALL RANGE (FIMAX, R, 2) 

404 R=R-DX( 1 1 ) 

400 CALL RANGE (R, FIMAX, 3) 

FO ( 1 , I ) =FTMAX 
GO TO 20 

12 IF (PROPNO(I) .NE. PROPNO ( 1 + 1 ) ) CALL PROPTY (2,1) 

FO(l,I)*FI(KFn 
DO 407 N*2,XEI 

I F ( PROPNO ( I ) .NE. PROPNO (I + D) CALL RANGE ( El ( N > ,DI ST ( N ) , 2 ) 
DIST(N)*DIST(N)+DX( I ) 

407 CALL RANGF(DIST(N) ,FI (N) ,3) 

GO TO (230,20) ,X0SW15 

C.....DFO-SMALL DFLTA F AT INTFRNAL PRINT POUNDS 
230 DEO * F0(1,I)/FNDE 
DO 9997 LL=1 ,NDF 
9997 E0(LL+1,I ) =E0 ( LL , I ) -DEO 
XOUNT ( I ) * NDE+1 
<M1 * NDE 
GO TO 30 

20 XONTUR*XNTR(XLTM) 

DO 202 L = l, XONTUR 

I F ( EO ( 1 , I ) .LE. E SPEC (L, KLIM) *1,0001) GO TO 204 

202 CONTINUE 
204 XOUNTI I >*L 

KM1=K0UNT( I )-l 
DO 208 JJ=1,KM1 
I ND= FOUNT ( I ) - JJ 
208 E0( JJ+1,I )=ESPFC( IND,KLIM) 

IF (J . NE. 1) GO TO 30 

XEI = 1 0670 

XL I M*2 

30 IF ( XM1 ,LE. 1 ) GO TO 234 

C CALCULATE INCIDFNT ENERGIES FROM ASSUMED EXIT ENERGIES 

DO 29 N=2,XM1 

XF I =XFT+1 0710 

CALL RANGE ( FO (N , I ) ,° ,2 ) 0720 

DIST (XEl ) =R+DX ( I ) 

29 CALL RANGE(DIST(KEI) ,FI (XFI ) ,3) 
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234 KEI*KET+1 

DI ST ( KF I ) =DX ( I) 

10 CALL RANGF(DIST(KET) »FI ( KFT ) » 3 ) 

GO TO (35,37) ,K0SW15 

35 DEO=EI (<FI)/FNDF 
DO 36 LL=1,NDE 
KEI = KE I +1 

36 FI (KEI )“EI (KEI-1 )-DFO 
El ( KEI 1*0,0 

GO TO 45 

37 KONTUR=KNTR(KLIM) 

DO 38 L*1»K0NTUR 

IF(FI(KEI) .LF. FS Dc T (L »KLIM) *1.0001 ) r,0 TO 39 

38 CONTINUE 

39 LL=L 
KM1 =LL— 1 

DO 40 L*1»KM1 

KEI*KEI+1 

ind*ll-l 

40 EHKFI )=ESPFC( IND.KLIM) 

45 DO 500 L=1,KEI 

IF (FT r L ) .LF. RNDLOW) GO TO 502 
500 CONTINUE 
GO TO 503 

502 KEI =L 

El (KEI )*BNDLOW 

503 READ (5*41) MOVF.EONO, TITLE 

41 FORMAT (2I3*11A6) 

MAx=KEI-l 

55 GO TO (62*64) *MOVE 

62 CALL FLUXEQ ( El ( 1 ) , OP ( 1 ) , EGNO ) 

64 DO 65 J=2,KEI 

DEI ( J-l ) = FI ( J-l )— FI ( J) 

C EIBAR ( J (-AVFRAGF ENFRGY 

EIBAR (J-l) = (PI ( J-l )+FI ( J) ) /2.0 
GO TO (142*145) ,MOVF 
142 CALL FLUXFQ ( FI ( J ) .OP ( J ) .FGNO ) 

C OPPRMf J)— INTEGRAL SPECTRUM 

OPPRM(J-l) * 0R(J)-0P(J-1) 

GO Tn 65 

145 CALL FLUXFQ ( EIBAR ( J-l ), OP ( J-l ) ,FQNO ) 

OPPRM( J-l )=0P( J-l )*DEI (J-l) 

65 CONTINUE 

GO TO (80,90) ,Knsvn 
90 WRITF (6,92) TITLE 
92 FDRMAT(1HL8X,11A6) 

L * 1 

GO TO (414*430) ,K0SW19 
414 DO 151 J=1 ,NOX 
I = N0XP1-J 
LIMIT = KOUNT ( I ) -1 
DO 425 M=l, LIMIT 

IF ( MOD ( L , 55 ) .NF. 1) GO TO (420,472) ,MOVF 
410 GO TO ( 412,415 ) *MOVF 

412 WRITF (6,413) 

413 F0RMAT(lH117X .8HENERGY , E8X , 17HN ( GREAT FR THAN F ) 8X .7HDFLTA F14X.6HE 
1 AVG.9X»15H(DN/DE)*DELTA E/2 IX , 3HMEV1 OX » 13HPR0T0NS/CM**214X .3HMEV1 
27X,3HMEV10X,13HPR0T0NS/CM**2) 

I F ( K0SW13 ,FQ. 1 )WRITF(6, 13413) 

13413 FORMAT ( 1H+ , 2 ( 46 X , 4H-S PC 1 OX ) ) 

GO TO 420 
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415 WRITE (6,416) 

416 FORMAT ( 1H117X * 8HENERGY »E1 3X * 7HDELT A E14X,6HF AVG. 10X , 1 3HDN/DF ( E AV 
1G. )6X,15H(DN/DE)*DPLTA F/12X,3(9X,3HMfV8X) , 17HPR0t0NS/CM**2"MEV5X , 
213HPROTDNS/CM**2 ) 

IF (KOSW13 ,FO. 1) WRITEt 6,16416 ) 

16416 FORMAT ( 1H+88X »4H— SFC14X »4H— SFC) 

422 WRITF (6,421) L ,E I ( L ) ,DFI ( L 5 , FIBAR ( L ) , OP ( L ) ,OPPRM ( L ) 

GO Tn 423 

420 WRITE (6,421) L , E I ( Lj ,0P ( L ) ,DFI ( L > ,ETBAR ( L ) ,OPPRM ( L ) 

421 FORMAT( I5,3X,1P5E20,5) 

423 IF ( L .GE. MAX) GO TO 443 
425 L = L+l 

151 WRITF (6,427) 

427 FORMAT ( 1H ) 

430 00 435 M*L,MAX 

IF ( MOD ( M , 55 ) ,NF. 1) GO TO t 439 ,441 ) , MOVP 

432 GO TO (453,436) , MOVE 

433 WRITE (6,413) 

IF ( K0SW13 .FQ. 1) WRITF(6, 13413 ) 

GO TO 439 
436 WRITE (6,416) 

IF (K0SW13 .FQ. 1) WRTTF(6* 16416) 

441 WRITE (6,155) M , El ( M ) ,OEI f M ) ,FT BAR ( M ) ,0P ( M ) ,OPPRM ( M ) 

GO TO 435 

439 WRITE (6,155) M ,E I ( M ) ,0P ( M ) ,DFI ( M) ,F IBAR ( M ) ,OPPRM ( M ) 

155 FORMAT ( I5,3X,1P5F20.5) 

435 CONTINUE 

443 GO TO (444,446) , MOVE 

444 WRITE (6,155) KFI,FI(KFT) ,OP(KFI ) 

GO TO 80 

446 WRITE (6,155) KFI »FT (KFT ) 

80 RETURN 
END 
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APPENDIX B 


SIBFTC TAPFIX 

DIMfNfi ON fNfRgy ( ion ) , R aNgf ; inn i , eNRgPR ( inn ,5 ) , ePRPR ( ion ,5 ) , 

A FPRNIK ion, 5 , .fnRRNIJ ( 1 n n , 5 ) , FNIJPR ( l ft 0 » 5 ) *FNI!NtJ(10fl»5) , 

P FUFRPR ( inn , «5 j , YPRCP ( 100 , 5 ) , yprc.N ( inn , e ) , vpppN ( inn , s ) , 

C FNFRNIM 100 , <5 ) ,YNUCP ( ion ,F) ,YNnCN( inn,*)) ,yN(ien (ino,5 ) # 

n FROMP ( inn ) ,XS DC> ( inn ! ,rpnMN ( inn ) ,ynMM( 1 no ) 

DIMENSION TABLFt 100*6) ,LFNGTH (4 ) ,RBENRG( inn ) ,RRF ( 100 ) ,K1NRG ( 100) , 

A <1 ( ion ) .K2NRG ( inn) ,K2 ( 100) , KNRG (100),KK(ino) , FXFNRG ( 100 ) 

B ,DFDX ( 1 O n ) 

EQUIVALENCE ( TABLF( 1 ,1 ) .K1NRG ( 1 ) ) , (TARLE( 1 ,2) ,Kl ( 1 ) ) , ( TABLE( 1 ,3) , 
AK2NRG ( 1 ) ) . (TABLES 1, A) .<2(1)5, (TABLFt 1 .5) .KNRG ( 1 ) ) , ( TARLF ( 1 ,6) .<<( 1 
B) ) 

. EQUIVALENCE ( LENGTH ( 1 ) ,L8 ) , (LENGTH ( 2 ) * L 9 ) , ( L p NO T H ( 3 ) ,L10 ) , ( LENGTH 
A (4) ,L11 ) 

REAL KlNRp,Kl,K2NRG,K?,KNR0»KK 
REWIND 3 

READ (5*4) N0MA T 1 .N0MAT2 
4 FORMAT ( 214) 

WRITF (3) NOMATi 
DO 2 N=l, NOMATi 

READ (5.6) MATNO.NOFCnM, GMW T ,L1 ,L2.L3,L4,L5,L6*L7 
6 F0RMA t ( 2 T 4 , F13, 6, 7T4) 

C RANGF-RNFRGy. 

READ (5,8) ( FNFRGY ( J ) ,RANG P '(J),J = 1,L1) 

8 FORMAT ( 8F9. 3 ) 

C ENERGY OF CASCADE PARTICLES, 

C A.) PROTONS BOMBARDING. 

DO 13 <=l,NOFCnM 

13 READ (5*14) (FNRGDR(J,<),epRPR{J,K),eprnu(J,<),J=1,L2) 

1 4 FORMAT ( 9F« . 2 ) 

C B.) NFUTRONS BOMBARDING. 

DO 15 K = 1 , NO FCOM 

15 READ (5*14) (FNRGNU(J.K) ,FNUPR( J,K) ,ENUNU( J,K) ,J=1,L3) 

C EMITTFO YIELDS. 

C A.) PROTONS BOMBARDING. 

DO 9 <=l,NnFC n M 

9 READ (5*10) (FNFRPRt J.K) ,YPRC D ( J,<) ,YPRCN( J.K) ,YPREN( J,K) ,J=i,L4) 

in ebrmat ( 8E9.0 ) 

C B.) NFUTRONS BOMBARDING. 

DO 11 < = 1 , NO FCOM 

11 Read ( 5 * i o > (enfRnu(J,o ,ynijcp( J.kj .ynucni j.o ,ynufN( j,k> ,j=i,L5) 

g X-SFCT I ONS . 

C A.) PROTON. 

RFAD (5.12) (FROMP(J) ,XSPR(J) ,J=1 »L6) 

12 FORMATf 10F7.0) 

G B.) NEUTRON. 

RFAD (5.12) (FPOMN(J) ,XSNU( J) ,J=1 ,L7) 

2 WRITE (3) MATNO.NOFCOM, GMWT ,L1»L2.L3,L4.L5,L6,L7, (FNERGY(J) .RANGF 
1 ( J ) , J= 1 . L 1 ) , ( ( FNRGPR (J.K) ,EPRPR ( J ,<) , fPRNU (J.K) ,J = 1,L2).< = 1 .NOFCOM 
2) . ( (ENRGNU(J.K) ,FNuPR(J,K) ,ENUNU(J,K) ,J = 1,L3) ,K = l .NOFCOM) ,( (FNFRPR 
3 ( J.K) , YPRCP ( J.K) , Y°RCN (J.K! , Y°REN (J.K) ,J=1 ,L4 ) ,K=1 .NOFCOM) , ( ( ENERN 
4U (J.K) .YNUCP ( J.K ) , YNUCNI J.K 5 , YNU EN ( J , K ) , J = 1 , L 5 ) ,K = 1 .NOFCOM) , ( FROMP 
s(J) »XSPR(J),J = i,L6 ).( erOMM ( J ) ,XSNIJ ( J ) , J= 1 ,L7 ) 

READ (5.18) L8.L9.L10.L11 

REAP (5.20) (RBFNRGf J) ,RRF( J) ,J=1 ,L8) 

READ (5.2 n ) (KlNRG(J) ,K1(J) »J=1 ,L9) 

RFAD (5.20) (K2NRG( J) ,K2( J) . J=1 ,L10) 

RFAD (5.20) (KNRG(J) , KK < J ) , J= 1 , L 1 1 ) 

DO 22 L = 1 » ’ 
m A x = L + L 

LIMIT=LENGTH(L+1 ) 
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DO 22 J=1,LI W IT 

22 TABLE ( J*MAX ) = TAPL r (J* w AX ) /360O . 

WRIT (3) LR*L9*L1 a ,L 11 , (RBENRE ( j) ,ppr { J ) » J= 1 ,L R ) , (KlNRr-(J) ,K1 (J) , 
AJ=1 *L9 ) * (K2NRO ( J) »K?( J) .J = l*LlO) » ( K.NRO { J ) »<K( J) *J = 1 *Lll ) 

V/R I TF (3) N'nVAT 
do 116 N=l,MOMAT2 
READ ( 3 * 1 R ) M AT"" 1 2 * L 1 2 
IP FORMAT (414) 

C DE/DX 

READ (3*20) (EXPNRr-(J) ,DFDX(J) ,J = 1,L12) 

2 A R0RMAT(F7. A ,F9.3*E7. A ,^9.3,P7. A ,F9.3, c 7, A ,c9.3) 

116 WRIT (3) MaTND?*L12» TXfNRC,! J) ,DPDX( J) »J=1*L12) 

END FILE 3 
REWIND 3 

CTflD 

END 
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APPENDIX C 


PROGRAM RUNNING TIME 


An order of magnitude estimate of the running time of the program may- 
be obtained as follows. Count the number of 5x's, 6E's, number of angles 
used in calculating evaporation neutrons, and the number of print bounds, 
then compute the following. 


Running time in minutes = 0.5 + 6.2X10"^ 



Np(5E)N p (6x)N p (LS) 


+ 0.08 



Np(PB) 


P = the number of problems stacked in the data deck 

Np(6E) = number of proton 5E’s in problem p (This has to be estimated 
1 on the basis of previous runs.) 

Np(6x) = number of proton Sx's in problem p 

Np(LS) = number of angles in problem p 

Np(PB) = number of print bounds in problem p 

This estimate does not reflect the time required to print spectra at each 
print bound. Also this estimate is based on calculating all generations of 
secondaries and the use of 69 energy groups for the cascade neutrons. If 
these conditions are changed new coefficients for 6.2X10“^ and 0.08 can 
be calculated by running a few cases. 


p = 1 

N p (5E) = 140 
Np(Sx) = 40 
Np(LS) = 5 
N p (PB) = 3 


Comparison of Running Time for Sample Problem 
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Running time in minutes (from formula) =2.5 min. 
Actual time was =2.3 min. 
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Thickness scale 


Figure 1. - Energy groups at internal bounds. Vertical 
lines show print boundaries. Tick marks on vertical 
lines show energy intervals at the print boundaries. 
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Figure 2 . - Construction of print boundaries 
and secondary source layers. 



